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Abstract 

Many  important  real-world  problems  require  meshing,  that  is  the  approximation  of  a  given 
geometry  by  a  set  of  simpler  elements  such  as  triangles  or  quadrilaterals  in  two  dimen¬ 
sions,  and  tetrahedra  or  hexahedra  in  three  dimensions.  Applications  include  finite  element 
analysis  and  computer  graphics.  This  work  focuses  on  the  former. 

A  physically-based  model  of  interacting  “particles”  is  introduced  to  uniformly  spread  points 
over  a  2-dimensional  polygonal  domain.  The  set  of  points  is  triangulated  to  form  a  triangle 
mesh.  Delaunay  triangulation  is  used  because  it  guarantees  a  low  computational  cost  and 
reasonably  well-shaped  elements.  Several  particle  interaction  (repulsion  and  attraction) 
models  are  investigated  ranging  from  Gaussian  energy  potentials  to  Laplacian  smoothing. 
Particle  population  control  mechanisms  are  introduced  to  make  the  size  of  the  mesh  elements 
converge  to  the  desired  size. 

In  most  applications  spatial  mesh  adaptivity  is  desirable.  Triangles  should  not  only  adapt  in 
size  but  also  in  shape,  to  better  fit  the  function  to  approximate.  Computational  fluid  dyna¬ 
mics  simulations  typically  require  triangles  stretched  in  the  direction  of  the  flow.  A  metric 
tensor  is  introduced  to  quantify  the  stretching.  The  triangulation  procedure  is  changed  to 
generate  “Delaunay”  meshes  in  the  Riemannian  space  defined  by  the  metric. 

This  new  approach  to  mesh  generation  appears  quite  promising. 
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Chapter  1 


Introduction 


Many  important  real-world  problems  require  meshing,  that  is  the  approximation  of  a  given 
geometry  by  a  set  of  simpler  elements  such  as  triangles  or  quadrilaterals  in  two  dimen¬ 
sions,  and  tetrahedra  or  hexahedra  in  three  dimensions.  Applications  include  finite  element 
analysis  [1]  and  computer  graphics. 

Many  engineering  simulations  require  the  solution  of  partial  differential  or  integral  equ¬ 
ations.  Since  most  of  these  equations  cannot  be  solved  analytically,  approximations  must  be 
used.  If  the  domain  has  a  simple  shape,  one  can  use  finite  difference  methods  with  structu¬ 
red  grids.  For  more  complex  domains,  finite  element  methods  are  used  on  an  unstructured 
grid,  that  is  a  mesh. 

Any  mesh  generator  should  address  the  following  concerns: 

•  functionality.  Obviously  it  should  work. 

•  robustness.  It  should  work  all  the  time,  with  any  input. 

•  quality.  The  quality  of  the  resulting  mesh  should  be  good,  that  is  closely  match  the 
desires  of  the  user. 

•  speed.  The  generation  process  should  be  fast. 

•  minimal  user  interaction.  Everything  that  can  be  automated  should  be. 

•  controllability.  The  user  should  be  able  to  influence  the  result  in  predictable  ways. 

In  this  work  we  limit  ourselves  to  generate  triangular  meshes  inside  a  two  dimensional 
domain  bounded  by  a  polygon  which  may  contain  holes.  The  main  application  we  focus  on 
is  finite  element  analysis. 

A  physically-based  model  of  interacting  “particles”  is  introduced  to  uniformly  spread 
points  over  the  domain.  The  set  of  points  is  triangulated  to  form  a  triangle  mesh.  De¬ 
launay  triangulation  is  used  because  it  guarantees  a  low  computational  cost  and  reasonably 
well-shaped  elements.  Several  particle  interaction  (repulsion  and  attraction)  models  are 
investigated  ranging  from  Gaussian  energy  potentials  to  Laplacian  smoothing.  Particle  po¬ 
pulation  control  mechanisms  are  introduced  to  make  the  size  of  the  mesh  elements  converge 
to  the  desired  size. 

In  most  applications  spatial  mesh  adaptivity  is  desirable.  Triangles  should  not  only 
adapt  in  size  but  also  in  shape,  to  better  fit  the  function  to  be  approximated.  Computational 
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fluid  dynamics  simuiations  typicaiiy  require  triangies  stretched  in  the  direction  of  the  flow. 
A  metric  tensor  is  introduced  to  quantify  the  stretching.  The  trianguiation  procedure  is 
changed  to  generate  “Deiaunay”  meshes  in  the  Riemannian  space  defined  by  the  metric. 


1.1  Previous  Work 

There  are  severai  surveys  avaiiabie  on  mesh  generation  [2,  15]. 

Most  of  present  mesh  generation  aigorithms  are  structured  in  the  foiiowing  way.  First  a 
mesh  is  buiid  with  methods  such  as  advancing  front  [19,  21,  23],  quadtree  decomposition  [35], 
or  by  greedy  point  insertion  [3,  32].  The  quaiity  of  the  mesh  is  further  improved  with  the 
use  of  smoothing.  The  most  common  method  is  Lapiacian  smoothing  [8]. 

Advancing  front  methods  start  meshing  at  the  boundaries  of  the  domain.  A  iist  of  nodes 
to  be  expanded  (referred  to  as  the  front)  is  maintained.  At  each  iteration,  the  front  advances 
by  expanding  a  node,  and  inserting  a  new  node  at  the  desired  distance  from  the  front.  Badiy 
shaped  eiements  can  appear  in  the  middie  of  the  domain,  where  the  fronts  coiiide. 

Quadtree  methods  recursiveiy  spiit  a  square  surrounding  the  domain  into  four  smaiier 
squares,  untii  the  desired  size  is  reached.  To  obtain  a  trianguiar  mesh,  the  squares  are 
further  divided  into  triangies. 

Greedy  insertion  methods  iocaiize  pooriy  shaped,  or  pooriy  sized  eiements,  and  spiit 
them  by  inserting  a  new  node  on  an  edge  [3] ,  at  the  center  of  a  triangie  [32] ,  or  at  the  center 
of  the  circie  circumscribing  a  triangie  [25].  Aithough  the  Deiaunay  criterion  is  the  most 
common,  other  triangie  quaiity  criteria  have  been  used  [31]  to  determine  the  topoiogy  of 
the  mesh. 

Lapiacian  smoothing  consists  of  moving  each  vertex  to  the  centroid  of  its  neighbors. 
This  operation  must  generaiiy  be  repeated  severai  times  for  each  node  before  the  quaiity 
of  the  mesh  is  improved.  Aithough  Lapiacian  smoothing  generaiiy  improves  the  shape  of 
eiements,  it  is  not  guaranteed  to  do  so,  especiaiiy  if  the  domain  is  concave.  Furthermore  the 
degrees  of  the  nodes,  which  are  often  the  reason  of  pooriy  shaped  eiements,  are  not  affected. 
Reiaxations  methods  that  combine  Lapiacian  smoothing  with  iocai  topoiogicai  optimization 
have  been  proposed  to  remedy  to  this  probiem  [10,  11]. 

Anisotropic  mesh  generation  has  not  benefited  from  an  extensive  iiterature.  One  of  the 
major  fleids  of  appiication  for  anisotropic  meshes  is  computationai  fluid  dynamics  (CFD), 
where  stretched  triangies  oriented  in  the  direction  of  the  flow  are  desirabie.  Mavripiis  [22] 
proposed  to  stretch  the  piane  by  iifting  it  on  a  surface  in  three  dimensions.  The  deformation 
of  space  is  represented  at  each  point  by  two  vaiues:  an  angie  giving  the  direction  of  the 
stretching,  and  a  vaiue  iarger  than  1  quantifying  the  stretching.  Before  that,  Peraire  [23] 
was  using  a  simiiar  representation  to  quantify  the  desired  eiement  size  as  a  function  of  its 
position  and  orientation. 

Later  a  metric  tensor  was  introduced  [3,31].  The  tensor  representation  has  the  advantage 
to  be  directiy  reiated  to  the  Hessian  of  the  function  (speed,  pressure,  etc.)  to  be  estimated  [5]. 
It  is  thus  a  more  natural  representation  to  create  adapted  meshes.  Quite  impressive  results 
have  been  produced  by  Castro-Diaz,  Hecht,  and  Mohammadi  [3].  In  their  method,  which  is 
of  the  greedy  insertion  type,  edges  that  are  too  long  are  split,  and  edges  that  are  too  short 
are  collapsed.  When  the  desired  number  of  elements  is  reached,  the  mesh  is  further  relaxed 
to  improve  its  quality. 

Mesh  generation  has  also  benefited  from  other  kinds  of  approaches.  The  idea  of  using 
physically-based  simulations  for  mesh  generation  has  been  investigated  by  Shimada  [26]  and 


his  bubble  packing  method.  The  bubble  interaction  (attraction/repulsion)  model,  which 
is  inspired  by  the  Lennard-Jones  interaction  model  from  molecular  chemistry,  generates 
triangulations  that  imitate  Nature  in  her  way  of  producing  regular  arrangements  of  points, 
such  as  in  crystals.  The  main  advantages  of  this  method  are  good  point  placement,  and 
intrinsic  remeshing  capabilities. 

Physically-based  models  have  also  been  used  in  computer  graphics  for  sampling  surfa¬ 
ces  [29,  30,  33,  34].  In  these  models,  particles  spread  over  complex  surfaces  to  form  uniform 
sampling  patterns.  Szeliski  [29]  used  attracting  and  repelling,  oriented  particles  to  intera¬ 
ctively  sculpt  surfaces.  Turk  [30]  considered  resampling  polygonal  surfaces  using  repelling 
particles.  Witkin  and  Heckbert  [34]  used  repelling  particles  to  sample  implicit  surfaces. 
Although  many  have  thought  of  introducing  anisotropy,  and  defining  the  sampling  density 
based  on  surface  curvature,  few  have  done  it  [30].  Witkin  and  Heckbert  [34]  have  noticed 
that  such  schemes  can  produce  very  regular  patterns  of  points.  Their  work  has  been  the 
starting  point  of  the  work  presented  in  this  document. 

1.2  Approach  Overview 

In  this  work,  we  limit  ourselves  to  generate  triangular  meshes  over  polygonal  domains  in 
the  plane.  The  domain  can  contain  holes,  and  constraints  such  as  line  segments  and  points. 
Constrained  points  define  nodes  that  should  appear  in  the  mesh,  and  constrained  line  se¬ 
gments  should  not  be  crossed  by  any  edge  in  the  mesh.  The  domain  can  be  represented  by 
a  planar  straight  line  graph  (PSLG),  which  is  a  set  of  points  and  non-crossing  edges.  An 
example  of  such  a  graph  is  given  in  Figure  1.1. 


Figure  1.1:  A  planar  straight  line  graph 

Given  such  a  domain,  points  are  added  inside  it  and  triangulated  to  form  a  mesh  such 
that  the  length  of  every  edge  matches  as  closely  as  possible  a  feature  size  function^ .  This 
function  is  given  in  input  to  the  mesher,  and  is  defined  in  terms  of  the  position  of  the  vertices 
of  an  edge.  Different  classes  of  such  functions  generate  different  kinds  of  meshes  such  as 
(see  Figure  1.2): 

constrained  where  no  feature  size  function  is  specified 
^We  also  refer  to  the  feature  size  function  as  the  desired  edge  length 
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CDT 


Uniform 


Isotropic 


Anisotropic 


Figure  1.2:  Different  kinds  of  meshes 


uniform  where  the  feature  size  function  is  a  constant 

isotropic^  where  the  desired  edge  length  depends  on  the  position  of  the  edge 

anisotropic  where  the  desired  edge  length  depends  on  the  position  and  the  orientation  of 
the  edge 

To  distribute  the  nodes  inside  the  domain,  a  physically-based  model  of  interacting  “par¬ 
ticles”  is  introduced.  The  functionality  of  the  model  can  be  expressed  as: 

1.  Input:  polygonal  domain  and  feature  size  function 

2.  build  an  initial  triangulation  of  the  domain  (constrained  Delaunay),  and  create  a 
particle  on  each  vertex  of  the  domain 

3.  create,  move,  and  annihilate  particles  inside  the  domain  until  equilibrium  is  reached. 
The  triangulation  of  the  set  of  particles  is  maintained  Delaunay  at  all  times,  that  is 

^  These  kind  of  meshes  are  also  referred  to  as  graded  meshes 
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the  topology  of  the  mesh  is  locally  optimized  after  each  particle  movement,  creation, 
and  annihilation. 

4.  Output:  a  nice  mesh 

To  illustrate  this  algorithm,  consider  the  domain  depicted  in  Figure  1.3a.  First  its 
constrained  Delaunay  triangulation  is  built  (Figure  1.3b). 

Then  the  physically-based  process  is  started.  At  each  step,  a  particle  is  randomly  picked, 
and  its  position  updated  according  to  the  positions  of  its  neighbors.  Then  the  local  particle 
density  is  estimated,  and  a  particle  is  created/annihilated  if  the  density  is  too  low/high. 
Figure  1.3c-f  shows  the  evolution  of  the  mesh.  In  a  first  phase  (approximatively  first  1000 
steps  in  this  case) ,  particles  are  created  until  the  population  reaches  the  desired  level  which 
depends  on  the  feature  size  function.  The  growth  is  regulated  by  an  adaptive  population 
control  scheme.  In  a  second  phase,  the  mesh  is  regularized,  that  is  node  placement  is 
improved.  For  this  example,  after  8410  steps,  equilibrium  is  reached  and  the  algorithm 
halts. 

During  the  whole  simulation  the  triangulation  is  maintained  Delaunay  using  procedures 
for  point  insertion,  motion  and  removal. 

1.2.1  Document  Outline 

The  outline  of  this  document  is  as  follows.  Chapter  2  introduces  Delaunay  triangulations 
and  related  algorithms.  Chapter  3  describes  the  physically-based  model  of  interacting  par¬ 
ticles.  Chapter  4  is  an  interlude  in  which  the  use  of  Java  for  programming  an  interactive 
Delaunay  triangulator  on  the  web  is  presented.  Chapter  5  generalizes  the  model  described 
in  Chapter  3  to  the  anisotropic  case,  and  presents  some  results.  Finally  conclusions  are 
drawn  in  Chapter  6. 
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a)  Input  PSLG  b)  Initial  triangulation  (CDT) 


c)  After  30  iterations  (64  particles)  d)  After  300  iterations  (226  particles) 


)  After  3000  iterations  (295  particles)  f)  After  8410  iterations  (295  particles) 


Figure  1.3:  Mesh  Evolution 


Chapter  2 

Delaunay  Triangulations 


The  Voronoi  diagram  (VD)  of  a  set  5  =  {si,  S2,  •  •  • ,  Sn)  of  points  in  the  plane,  called  sites, 
is  a  partition  of  the  plane  into  n  convex  regions,  one  per  site.  Each  Voronoi  cell  V,  contains 
all  the  points  in  the  plane  closer  to  s,  than  to  any  other  site.  The  planar  dual  of  the  Voronoi 
diagram,  obtained  by  adding  a  line  segment  between  each  pair  of  sites  of  S  whose  Voronoi 
regions  share  an  edge,  is  called  the  Delaunay  triangulation  (DT). 

More  practical  definitions  of  a  Delaunay  triangulation  are  (all  of  the  following  statements 
are  equivalent): 

•  if  a  and  b  are  input  points,  the  DT  contains  the  edge  {a,  b}  if  and  only  if  there  is 
a  circle  through  a  and  b  that  intersects  no  other  input  points  and  contains  no  input 
points  in  its  interior 

•  the  circumscribing  circle  of  each  triangle  contains  no  input  points  in  its  interior. 

There  is  also  a  nice  relationship  between  Delaunay  triangulations  and  3-dimensional 
convex  hulls.  Lift  each  point  of  the  input  to  a  paraboloid  by  mapping  the  point  {x,y)  to 
{x,y,x^  +y‘^y  ■  It  can  be  proved  [7]  that  the  DT  of  the  input  points  is  the  projection  of  the 
lower  convex  hull  onto  the  a;y-plane. 

The  Delaunay  triangulation  features  other  interesting  properties.  Indeed  it  maximizes 
of  the  minimum  angle,  minimizes  the  maximum  circumcircle  as  many  other  measures. 

The  outline  of  the  chapter  is  as  follows.  Section  2.1  defines  data  structures  for  repre¬ 
senting  triangulations.  Associated  topological  and  geometrical  operators  are  defined  in  Se¬ 
ctions  2.2  and  2.3.  Section  2.4  presents  some  methods  for  building  Delaunay  triangulations. 
Section  2.5  addresses  the  problem  of  dynamically  maintaining  a  triangulation  (incremen¬ 
tal  site  insertion  and  removal).  Section  2.6  presents  constrained  Delaunay  triangulations 
(CDT). 


^This  property  also  holds  for  any  paraboloid  z  =  a-  ((x  —  a)^  +  {y  —  6)^)  where  a,  a,  and  b  are  constants. 
This  trivially  follows  from  the  fact  that  a  DT  is  invariant  to  translations  of  the  input  set  of  sites.  The  value 
of  a  doesn’t  matter  either  since  it  neither  affects  the  topology  of  the  convex  hull,  nor  its  projection  on  the 
plane 
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2.1  The  Quadedge  Data  Structure 

The  quadedge  data  structure  [13]  was  designed  for  representing  general  subdivisions  of  orien- 
table  manifolds.  It  simultaneously  represents  both  the  subdivision  and  its  dual.  Alternatives 
are  the  double-connected-edge-list  (DCEL)  [24]  and  the  winged-edge  data  structures,  which 
do  not  hold  the  dual.  The  quadedge  structure  has  been  preferred  because  an  implementation 
was  readily  available  [20]. 


e 

e.invRot 

e.Rot 

1 

e.Sym 

Figure  2.1:  The  four  directed  edges  of  a  quadedge 

Each  quadedge  record  holds  four  directed  edges  corresponding  to  a  single  undirected 
edge  in  the  subdivision  and  to  its  dual  edge.  Each  directed  edge  has  three  references:  Org, 
which  points  to  the  site  at  its  origin,  Rot,  which  points  to  its  dual  edge,  and  Onext,  which 
points  to  the  counterclockwise  next  edge  in  the  subdivision.  All  other  topological  operators 
(see  Figures  2.1  and  2.2)  can  be  defined  in  terms  of  these  primitives,  as  summarized  below. 


Figure  2.2:  Edge  navigation  operators 

Navigation  inside  a  quadedge  is  achieved  with  the  use  of  the  identity.  Rot,  Sym,  and 
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Rot  ^  operators,  all  of  which  can  be  expressed  in  terms  of  Rot  operations: 


e  =  e 

e.Rot  =  e.Rot 

e.Sym  =  e.Rot^ 

e.Rot~^  =  e.Rot^ 


(2.1) 


where  e.Rot  means  that  operator  Rot  is  applied  to  edge  e,  and  e.Rot"^  that  Rot  is  applied 
n  times. 

The  position  of  each  site  in  the  quadedge  can  be  retrieved  with  the  Org  operator: 


e.Org 
e.Left 
e.Dest 
e. Right 


e.Org 

e.Rot~^  .Org 

e.Sym.Org 

e.Rot.Org 


e.Org 

e.Rot^.Org 

e.Rot^.Org 

e.Rot.Org 


(2.2) 


where  the  middle  column  represents  the  definition  of  the  operator,  and  the  right  column  its 
expression  in  terms  of  primitives. 

Onext  allows  navigation  through  the  subdivision.  Movements  to  neighboring  edges  are 
represented  in  figure  2.2.  They  can  be  subdivided  into  two  classes,  namely  counter-clockwise 
ones: 


e.Onext 

e.Lnext 

e.Rnext 

e.Dnext 

and  clockwise  ones: 

e.Oprev 

e.Lprev 

e.Rprev 

e.Dprev 


e.Onext 

e.Rot~^  .Onext. Rot 
e.Rot. Onext. Rot~^ 
e.Sym. Onext. Sym 


e.Rot. Onext.  Rot 
e.Onext. Sym 
e.Sym.Onext 
e.Rot~^  .Onext. Rot~^ 


e.Onext 

e.Rot^  .Onext. Rot 
e.Rot. Onext.  Rot^ 
e.Rot^  .Onext. Rot^ 


e.Rot. Onext.  Rot 
e.Onext. Rot^ 
e.Rot^  .Onext 
e.Rot^  .Onext.  Rot^ 


(2.3) 


(2.4) 


When  assuming  that  the  subdivision  is  a  triangulation,  it  is  possible  to  simplify  some  of  the 
operations: 

e.Dprev  =  e. Onext. Rot^ .Onext 

e.Dnext  =  e.Rot.Onext^ .Rot  '  '  ' 

The  next  section  presents  higher  level  operators. 


2.2  Topological  Operators 

Starting  with  the  operator  set  proposed  by  Guibas  and  Stolfi  [13],  we  have  modified  it  to 
make  it  more  symmetrical,  and  also  to  eliminate  superfluous  edge  navigation.  The  parame¬ 
ters  of  the  Ccmnect  operator  have  been  changed  and  its  inverse  operator  Disccmnect  is  intro¬ 
duced.  The  latter  accommodates  some  operations  that  were  previously  part  of  DeleteEdge. 
As  a  consequence  DeleteEdge  has  also  been  changed.  Each  operator  is  described  in  the 
next  paragraphs. 


15 


2.2.1  MakeEdge 

MakeEdge  creates  a  new  quadedge  and  initializes  all  of  its  four  directed  edges.  It  takes  no 
argument  and  returns  the  first  edge  of  the  quadedge.  The  quadedge  is  a  subdivision  by 
itself,  namely  the  one  of  a  sphere  [13].  The  details  of  the  procedure  are  described  below. 
The  Create  operator  creates  a  new  directed  edge  and  initializes  all  its  pointers  to  nil. 

begin  Edge. MakeEdge 

this.Create,  62. Create,  ez.Create,  64. Create 

this. Onext  this, e2.0next  ei,  es.Onext  e^,  e^.Onext e2 
this. Rot  e2,  e2.i?ot  t— es,  e^.Rot  e^,  e 4,. Rot this 

end 

2.2.2  Splice 

Splice  is  the  basic  operator  used  to  attach  and  detach  edges  from  each  other  (see  figure  2.3). 
It  takes  an  edge  b  as  parameter  and  returns  no  value.  It  is  its  own  inverse.  Its  description 
is  given  below. 


Figure  2.3:  The  effect  of  the  Splice  operator 

begin  Edge.Splice{h) 

Cl  t—  this.Onext,  €2  •<—  h.Onext 

a  t—  Cl. Rot,  l3  t—  62. Rot 

63  t—  a.Onext,  64  t—  jS.Onext 

this .Onext  t—  62,  h.Onext  t—  ei 

a.Onext  t—  64,  jS.Onext  t—  63 

end 

2.2.3  Connect 

The  Connect  operator  connects  the  two  vertices  of  an  edge  e  to  respectively  two  edges  a 
and  b.  It  takes  as  argument  the  two  edges  a  and  b,  and  returns  no  value.  The  Ccmnect 
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operator  is  basically  a  succession  of  two  Splice  operations  followed  by  an  update  of  the  Org 
fields.  A  description  is  given  below. 

begin  Edge.Ccmnect{a,b) 
this.Splice(a) 
this.Org  t—  a.Org 
this.Sym.Splice{h) 
this.Dest  t—  b.Org 
end 

This  new  definition  of  the  Connect  operator  changes  in  two  ways  from  its  original  [13] 
form:  it  does  not  create  a  new  edge  any  more,  and  the  arguments  represent  different  edges 
so  that  they  can  directly  be  applied  to  both  of  the  Splice  operations. 

2.2.4  Disconnect 

The  Disconnect  operator  is  the  inverse  of  the  Connect  operator.  Thus^ 

e.Connect(a,b).Disconnect(a,b)  =  e 

Its  description  is  given  below.  Since  the  Splice  operator  is  its  own  inverse  the  Disconnect 
operator  is  quite  the  same  as  Connect.  The  difference  lies  in  the  update  of  the  Org  fields 
(which  are  reset  to  nil  by  Disconnect). 

begin  Edge.Disconnect(a  t—  this.Oprev,  b  t—  this.Lnext) 
this.Splice(a) 
this.Org  t—  nil 
this.Sym.Splice(b) 
this.Dest  t—  nil 
end 


2.2.5  Swap 

The  Swap  operator  is  at  the  core  of  the  edge  swapping  algorithm  described  in  section  2.4.1. 
To  swap  an  edge  is  to  replace  it  by  the  other  diagonal  of  the  quadrilateral  in  which  it  is 
inscribed.  Swap  is  its  own  inverse  and  thus  e. Swap. Swap  =  e. 

begin  Edge. Swap 
a  t—  this.Lprev 
b  t—  this.Sym.Lprev 
c  t—  a.Lprev 
d  t—  b.Lprev 
this .Disconnect(d,  c) 
this.Ccninect{b,  a) 
end 


^The  relation  only  holds  if  e.Org  =  e.Dest  =  nil 
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2.2.6  DeleteEdge 

The  DeleteEdge  operator  undoes  everything  the  MakeEdge  operators  does,  as  described 
below. 


begin  Edge. DeleteEdge 

62  t—  this. Rot,  63  t—  e2.Rot, 

t his. Rot  nil,  e2.Rot-^nil, 

this.Onext  nil,  02. 0 next  nil, 

this. Destroy ,  02. Destroy, 

end 


64  t—  e^.Rot 
e^.Rot  t—  nil, 
e^.Onext  t—  nil, 
ez. Destroy , 


04. Rot  t—  nil 
eA.Onext  t—  nil 
04. Destroy 


Putting  references  to  nil  and  then  Destroy-ing  edges  clearly  is  redundant.  The  reason 
both  appear  in  the  pseudocode  it  that  languages  such  as  C++  require  the  latter,  whereas 
languages  featuring  automatic  garbage  collection  require  the  former. 


2.3  Geometrical  Operators 

Geometrical  operators  are  defined  to  test  the  relative  positions  of  points  and  edges.  There  is 
also  an  InCircle  operator  which  tests  whether  a  point  lies  inside  the  circle  defined  by  three 
other  points,  and  which  is  one  of  the  most  important  operators  for  constructing  Delaunay 
triangulations. 

2.3.1  OnRight,  OnLeft  and  OnEdge 

The  OnRight,  OnLeft  and  OnEdge  operators  tell  if  a  given  point  is  respectively  to  the 
right  of,  to  the  left  of,  or  on  an  edge.  The  decision  is  based  on  the  sign  of  the  signed  area 
of  the  triangle  defined  by  the  two  vertices  of  the  edge  and  the  point  to  test. 

begin  Edge.OnRight(p) 

return  (p  —  this.Org)  x  {this.Dest  —  this.Org)  >  0 

end 

begin  Edge.OnLeft{p) 

return  (p  —  this.Org)  x  {this.Dest  —  this.Org)  <  0 

end 

begin  Edge.OnEdge{p) 

return  (p  —  this.Org)  x  {this.Dest  —  this.Org)  =  0 

end 

where  x  is  the  cross  product  between  two  vectors. 

2.3.2  InCircle 

Given  an  edge  e  and  two  points  i  and  j,  the  InCircle  operator  returns  true  if  and  only  if 
the  circle  defined  by  i  and  the  two  vertices  of  e  includes  j.  It  tells  whether  an  edge  should 
be  swapped  or  not.  This  test  is  at  the  heart  of  most  Delaunay  triangulation  algorithms,  and 
is  floating  point  intensive.  In  practice  it  is  important  to  optimize  it.  Three  different  ways 
of  performing  the  test  are  presented,  and  their  speed  is  evaluated. 
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Circle  Test 


The  naive  way  is  to  compute  the  position  of  the  center  and  the  radius  of  the  circle,  and  then 
check  if  the  distance  from  its  center  to  j  is  less  than  its  radius.  The  center  of  the  circle  lies 
at  the  intersection  of  the  two  lines  that  respectively  bisect  e.Org  and  i,  and  e.Dest  and  i. 
These  two  lines  can  be  expressed  as 


rii  ■  X  —  Cl  =  0 
712  •  X  -  C2  =  0 

where  rii  =  e.Org  —  iis  a  normal  to  the  first  line  and  ci  =  ^(e.Org  +  i)  •  rii  is  a  constant. 
712  and  C2  are  similarly  defined  for  the  second  line.  The  two  equations  can  be  combined  into 
a  single  matrix  equation  Nx  —  c  =  0.  The  solution  of  which  is  a;  =  N~^c,  where  x  is  the 
intersection  of  the  two  lines,  and  thus  the  center  of  the  circle. 

The  InCircle  test  then  summarizes  to 

\x-  3\  <\x-i\ 

To  compute  the  center  of  the  circle,  11  multiplications,  13  additions,  and  2  divisions 
are  required.  The  distance  test  requires  additional  6  additions  and  4  multiplications.  The 
total  cost  is  thus  15  multiplications,  19  additions,  and  2  divisions.  It  is  possible  to  avoid 
the  divisions  at  the  cost  of  four  extra  multiplications  (by  multiplying  everything  by  the 
determinant  of  N). 


Convex  Hull  Test 


A  better  way  might  be  to  consider  some  other  properties  of  Delaunay  triangulations.  The 
convex  hull  test  is  based  on  the  nice  relation  between  Delaunay  triangulations  and  3- 
dimensional  convex  hulls.  The  test  is  based  on  the  sign  of  the  signed  volume  of  the  tetra¬ 
hedron  defined  by  the  projection  of  the  four  points  on  a  paraboloid: 


Xi 

yi 

xj  +y‘f 

1 

X2 

X3 

yi 

yz 

xl  +yl 
xl  +y‘l 

1 

1 

>  0 

X4 

yi 

xl  +yl 

1 

=  t 

z.Dest,  {X3,y3) 

= 

i,  and 

(2.6) 


,^4)  =  j.  When  naively 
computed,  this  determinant  requires  48  multiplications  and  27  additions.  A  good  idea  is  to 
consider  that  the  last  column  is  filled  with  ones  to  reduce  the  cost  to  20  multiplications  and 
24  additions,  as  in  the  code  proposed  by  Lischinski  [20]. 

It  is  possible  to  improve  this  result  by  translating  all  the  points  by  a  same  amount,  so 
that  one  of  the  points  lies  at  the  origin  (0,  0): 


0 

0 

0 

X2 

xl  +yl 

X3 

yz 

xl  +yl 

X4 

yi 

xl  +yl 

19 

(2.7) 


where  Xi  =  Xi  —  xi,  yi  =  yi  —  yi  and  the  paraboloid  is  Zi  =  (a;,  —  a;i)^  +  (y,  —  yi)^.  The  test 
is  thus  equivalent  to: 


X2 

h 

xl  +yl 

X3 

+yl 

>  0 

(2.8) 

Xi 

n 

xl  +yl 

which  requires  only  15  multiplications  and  14  additions. 


Angle  test 

It  also  possible  to  define  a  test  based  on  angles^:  if  the  sum  of  angles  a  and  /3  (see  figure  2.4) 
is  smaller  than  180  degrees,  then  the  edge  should  be  swapped.  Another  way  to  put  it  is 
sin(Q:  +  P)  >0?  This  expression  can  be  simplified  in  the  following  way: 


sin(Q:  +  P) 


=  sin  a  cos  /3  +  sin  /3  cos  a 

a  X  b  c  •  d  c  X  d  a  •  b 

|a|  •  |6|  |c|  •  \d\  |c|  •  \d\  |a|  •  |6| 

oc  {a  X  b){c  ■  d)  +  {c  X  d){a  ■  b) 


(2.9) 


This  test  only  requires  10  multiplications  and  13  additions. 


Figure  2.4:  The  InCircle  test 


Summary 

From  the  several  ways  of  performing  the  InCircle  test  that  have  been  examined,  the  angle 
test  is  fastest,  followed  by  the  convex  hull  test,  and  the  circle  test.  The  operator  is  thus 
defined  using  the  angle  test: 

^Thanks  to  Marshall  Bern  for  pointing  it  out 
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begin  Edge.InCircle{p  •<—  e.Onext.Dest,q  <r-  e.Oprev.Dest) 
a  q  —  this.Org 

h  p  —  this.Org 

c  p  —  this.Dest 

d  q  —  this.Dest 

return  (a  x  b){c  ■  d)  +  {c  x  d){a  ■  b)  >0 

end 

2.4  Triangulation  Algorithms 

Several  methods  are  known  for  triangulating  the  convex  hull  of  a  set  of  points,  namely  the  flip 
algorithm,  the  incremental  algorithm,  the  divide-and-conquer  algorithm  and  the  sweepline 
algorithm.  The  first  two  are  the  topic  of  the  next  sections.  The  divide-and-conquer  and 
sweepline  algorithms  are  only  briefly  described  because  they  haven’t  been  used  in  the  frame 
of  this  work. 

In  the  divide-and-conquer  algorithm,  the  input  site  set  is  first  sorted  by  a;-coordinate 
and  split  vertically  in  two  subsets  of  equal  size.  The  Delaunay  triangulation  is  computed 
recursively  for  each  subset.  Subsets  are  then  merged  using  a  sweeping  circle  algorithm.  A 
complete  description  is  given  in  [13].  The  time  complexity  of  the  algorithm  is  O(nlogn). 

The  sweepline  algorithm  is  due  to  Fortune  [9].  In  a  comparison  with  divide-and-conquer, 
Leach  [18]  has  optimized  both  algorithms  and  measured  execution  speed.  According  to 
Leach,  the  sweepline  algorithm  is  the  fastest,  but  also  appears  to  be  the  least  robust. 

2.4.1  The  Edge  Swapping  Algorithm 

Given  an  initial  triangulation,  the  edge  swapping  algorithm  maintains  a  queue  of  edges  that 
might  fail  the  circumcircle  test.  An  edge  e  is  said  to  pass  the  circumcircle  test  if  the  circle 
through  {e.Org,e. Best, e.Onext. Best}  does  not  include  the  point  e.Oprev.Dest.  In  the 
initial  triangulation,  any  edge  might  fail  the  test,  so  the  queue  initially  contains  all  edges. 
Then,  the  first  edge  e  is  repeatedly  removed  from  the  queue.  If  e  does  pass  the  circumcircle 
test,  the  next  edge  is  dequeued.  But  if  e  does  fail  the  test,  it  is  swapped.  As  this  operation 
might  change  the  status  of  some  of  the  four  edges  of  the  quadrilateral  of  which  e  is  a  diagonal, 
those  edges  are  added  into  the  queue  if  not  already  there.  The  algorithm  stops  when  the 
queue  is  empty. 

The  algorithm  always  terminates  after  O(n^)  flips.  Combined  with  an  O(n^)  algorithm 
for  constructing  the  initial  triangulation,  the  Delaunay  triangulation  of  a  set  of  points  can 
be  found  in  O(n^)  time  with  this  method. 

2.4.2  The  Incremental  Algorithm 

The  incremental  algorithm  was  first  proposed  by  Lawson  [17].  It  starts  with  a  triangle 
large  enough  to  contain  all  the  points  of  the  input  set  (ideally  the  three  vertices  are  at 
infinity).  Points  are  added  into  the  triangulation  one  by  one,  maintaining  the  invariant  that 
the  triangulation  is  Delaunay.  First  the  triangle  T  containing  the  new  point  p  is  located. 
New  edges  are  created  to  connect  p  to  the  vertices  of  T.  The  edges  of  T  are  inspected  to 
verify  that  they  still  satisfy  the  circumscribing  circle  condition.  If  the  condition  is  satisfied 
the  edge  remains  unchanged.  If  it  is  violated  the  offending  edge  is  swapped.  In  this  case 
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two  more  edges  become  candidates  for  inspection.  The  process  continues  until  no  more 
candidates  remain. 

If  the  order  of  insertion  is  randomized,  and  with  use  of  appropriate  data  structures 
for  point-in-triangle  location,  it  can  be  shown  [14]  that  the  expected  running  time  of  the 
algorithm  is  O(nlogn).  However,  without  any  additional  data  structure,  jump-and-walk 
methods  can  be  used  to  reach  an  expected  time  of  [6]. 

2.5  Dynamic  Maintenance  of  a  Triangulation 

The  problem  of  dynamically  maintaining  the  triangulation  can  be  divided  into  two  sub¬ 
problems:  insertion  and  removal  of  a  site.  Moving  a  site  can  be  seen  as  a  removal  followed 
by  an  insertion. 

2.5.1  Site  Insertion 

The  procedure  for  inserting  a  point  i  runs  as  follows: 

1.  locate  the  triangle  T  in  which  i  lies 

2.  add  an  edge  from  i  to  each  vertex  of  T 

3.  put  the  three  edges  of  T  on  a  queue 

4.  while  the  queue  is  non-empty,  remove  the  first  edge  e  from  the  queue.  If  e  does  pass  the 
circumcircle  test,  the  next  edge  is  dequeued.  If  e  fails  the  test,  swap  it,  and  insert  in 
the  queue  the  two  edges  a  and  b  adjacent  to  e  which  do  not  contain  i  (see  Figure  2.5). 


e.Swap 


Figure  2.5:  Edge  swapping  during  point  insertion 

It  is  usually  assumed  that  points  are  in  general  position,  that  is  that  no  three  points  are 
collinear  or  four  points  cocircular.  In  practice  however  this  is  not  the  case.  As  a  result  the 
insertion  of  a  site  i  on  an  edge  e  (three  colinear  point)  has  to  be  treated  in  a  special  way. 

1.  remove  the  edge  e.  i  then  lies  inside  a  quadrilateral  Q 

2.  add  an  edge  from  the  i  to  each  vertex  of  Q  which  is  the  containing  quadrilateral 

3.  put  the  four  edges  of  Q  on  a  queue 

Step  4  is  identical  to  the  one  for  general  position  point  insertion.  A  more  detailed 
description  of  the  point  insertion  routine  is  given  in  figure  2.6. 

The  number  of  edge  flips  is  equal  to  the  degree  of  i  after  its  insertion.  As  the  average 
degree  of  a  vertex  in  a  planar  graph  is  less  than  six,  the  average  cost  of  step  4  is  0(1).  Worst 
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begin  Insert(p) 
e  <r-  Locate(p) 

if  e.OnEdge{p)  then 
e  <r-  e.Oprev 
e.Onext.DeleteEdge 

base  <r-  MakeEdge 
base.Splice(e) 
base.Org  <r-  p 
base. Best  <r-  p 

b  <r-  base.Sym 

loop 

a  <r-  base.Lprev 
if  a  =  e  then 
exit  loop 

base  <r-  MakeEdge 
base.Ccmnect(a,  b) 

e  <r-  base 

do 

a  <r-  e.Onext 
b  <r-  a.Onext 

while  a.InCircle{b.Dest,p)  do 
a. Swap 
a  b 

b  <r-  a.Onext 
e  <r-  e.Dprev 
while  e  ^  base 
end 


Figure  2.6:  Point  insertion  procedure 


case  is  0{n)  and  happens,  for  example,  when  i  is  inserted  at  the  center  of  n  points  lying  on 
a  circle. 

The  point  location  cost  (step  1)  can  be  reduced  to  O(logn)  time  with  appropriate  0{n)- 
space  data  structures  [14].  Without  additional  data  structures  a  simple  walking  method 
can  be  used  to  achieve  an  expected  0(i/n)  time  performance.  A  fast  walking  method  is 
described  in  the  next  section. 

The  total  cost  of  an  insertion  is  thus  dominated  by  the  cost  of  the  point  location  proce¬ 
dure. 

2.5.2  Walking  Method  for  Point  Location 

Guibas  and  Stolfi  [13]  have  proposed  a  simple  walking  method  for  point  location,  which 
Lischinski  has  implemented  [20].  An  improved  version  where  redundant  tests  have  been 
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removed  is  described  by  a  finite  state  automaton  in  Figure  2.8.  At  any  time  the  point  to 
locate  is  on  the  left  of  the  current  edge.  Figure  2.7  shows  the  two  possible  movements  at 
each  step.  Performance  is  compared  in  Table  2.1  in  terms  of  function  calls.  It  turns  out 
that  the  new  algorithm  is  more  than  twice  as  fast  as  the  original  one. 


10000  locations 

100000  locations 

Procedure 

Guibas 

Bossen 

ratio 

Guibas 

Bossen 

ratio 

Edge. Righto  f 

207.8 

90.7 

2.29 

657 

277 

2.37 

Point.  = 

167.7 

78.1 

2.15 

531 

243 

2.18 

Edge.Onext 

126.1 

44.1 

2.86 

403 

137 

2.94 

Edge.Dprev 

80.1 

43.6 

1.84 

251 

137 

1.83 

Edge.Org 

83.9 

39.1 

2.15 

265 

122 

2.17 

Edge.Dest 

83.9 

39.1 

2.15 

265 

122 

2.17 

Edge.Sym 

0.52 

0.62 

0.84 

0.40 

0.62 

0.65 

TotaU 

374.5 

166.5 

2.25 

1183 

519 

2.27 

Table  2.1:  Guibas  versus  Bossen  point  location  performance  (listing  the  number  of  calls  to 
each  procedure,  per  insertion) 


Figure  2.7:  The  two  possible  walking  directions 


2.5.3  Site  Removal 

Unfortunately  removing  a  site  i  is  not  as  easy  as  inserting  one.  Although  it  can  be  simple 
when  using  appropriate  data  structures  and  removing  the  last  inserted  point  [16],  the  general 
case  is  more  difficult  to  handle.  The  procedure  description  is  still  very  short: 

1.  remove  all  the  edges  e  such  that  e.Org  =  i 

2.  delete  i 

3.  triangulate  the  simple  polygon  which  contained  i 
■^Sum  of  all  topological  operations  calls 
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begin  Locate(p) 

sO:  e  •<—  some  edge  /  such  that  ^f.RightOf{p) 

goto  if  e.Org  =  pV  e.Dest  =  p  then  s7  else  si 

si:  /  t—  e.Onext 

goto  if  f  .Righto f  ip)  then  s4  else  s6 
s2:  /  t—  e.Dprev 

goto  if  f .Righto f  ip)  then  s3  else  s5 
s3:  /  t—  e.Onext 

goto  if  f  .Righto f  ip)  then  s7  else  s6 
s4:  /  t—  e.Dprev 

goto  if  f  .Righto f  {p)  then  s7  else  s5 
s5:  e  t—  / 

goto  if  e.Org  =  p  then  s7  else  si 
s6:  e  t—  / 

goto  if  e.Dest  =  p  then  s7  else  s2 

s7:  retnrn  e 
end 


Figure  2.8:  Triangle  location  procedure 


The  cost  of  a  deletion  is  O(n^)  (see  next  section  on  simple  polygon  triangulation),  where 
n  is  the  degree  of  site  i.  As  the  degree  i  is  expected  to  be  less  than  six  on  average,  the 
asymptotic  time  is  not  relevant  in  practice. 

2.5.4  Triangulation  of  a  Simple  Polygon 

We  consider  here  the  case  of  a  hole  (simple  polygon)  inside  a  triangulation  that  needs  to 
be  retriangulated.  Triangulating  the  simple  polygon  is  not  as  easy  as  one  would  like  it 
to  be  since  the  polygon  cannot  be  assumed  to  be  convex.  A  simple  solution  is  to  build 
a  triangulation  of  the  polygon  by  successively  cutting  off  ears®,  and  then  run  the  edge 
swapping  algorithm. 

The  time  cost  of  the  retriangulation  is  O(n^),  where  n  is  the  number  of  sites  of  the  hole. 


2.6  Constrained  Delaunay  Triangulations 

A  constrained  Delaunay  triangulation  (CDT)  is  very  similar  to  a  Delaunay  triangulation. 
The  difference  is  that  some  edges  are  fixed  (constrained).  They  cannot  be  swapped  even  if 
they  do  not  pass  the  InCircle  test. 

®An  ear  is  a  triangle  which  shares  two  edges  with  the  polygon  to  be  triangulated 
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The  methods  for  generating  a  constrained  triangulation  are  quite  similar  to  the  ones  used 
for  unconstrained  triangulations.  The  edge  swapping  algorithm  for  instance  starts  with  a 
triangulation  where  all  constrained  edges  are  present,  then  runs  as  usual  but  a  constrained 
edge  is  never  swapped. 

There  are  also  0(n  log  n)  algorithms  for  generating  constrained  Delaunay  triangulation 
such  as  divide-and-conquer  [4],  but  they  are  usually  difficult  to  implement.  The  next  section 
presents  an  asymptotically  slower  incremental  algorithm  which  is  easier  to  implement. 

2.6.1  An  Incremental  CDT  Algorithm 

The  present  algorithm  for  building  constrained  Delaunay  triangulations  can  be  split  in  two 
phases.  In  the  first  one,  the  constraints  are  ignored,  and  a  triangulation  is  built  using  the 
incremental  algorithm  described  in  section  2.4.2.  In  the  second  phase,  the  constrained  edges 
are  inserted  one  by  one,  if  not  already  present.  First  all  the  edges  crossing  the  constrained 
edge  are  removed,  then  the  constrained  edge  is  inserted,  and  finally  the  simple  polygons 
on  each  side  of  the  constrained  edge  are  retriangulated.  The  procedure  for  removing  the 
crossing  edges  and  inserting  the  constrained  edge  is  described  in  Figure  2.9. 

begin  InsertConstraint(c) 

find  an  edge  e  such  that  e.Org  =  c.Org 
while  ^{e.OnLeft{c.Dest)  V  ^e.Onext.OnLeft{c.Dest)) 
e  t—  e.Onext 
eo  t-  e 
do 

/  e.Lnext 

while  c.OnLeft{f.Dest)  do 
f  .Remove 
f  t—  e.Lnext 
e  •<-  / 

while  e.Dest  ^  c.Dest 
e  t—  e.Lnext 
c.Splice{eo) 
c.Sym.Splice(e) 

retriangulate  the  hole  on  the  left  of  c 
retriangulate  the  hole  on  the  right  of  c 

end 


Figure  2.9:  The  InsertCcmstraint  procedure 
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Chapter  3 


A  World  of  Particles 


This  chapter  describes  a  system  of  interacting  particles  that  is  used  to  uniformly  spread 
points  over  a  domain.  At  any  time  the  set  of  points  defined  by  the  positions  of  the  particles 
is  triangulated  using  the  Delaunay  criterion.  The  goal  is  to  generate  a  particle  distribution 
such  that  every  edge  in  the  mesh  has  a  length  equal  to  a,  which  is  a  constant.  As  this 
optimality  criterion  is  usually  not  achievable  because  of  constraints  at  the  boundaries  of  the 
domain,  we  try  to  get  as  close  to  it  as  possible.  A  particle  interaction  scheme  is  shaped 
according  to  following  objectives: 

•  Simplicity.  The  system  should  be  governed  by  a  small  set  of  simple  rules. 

•  Speed.  The  computational  cost  for  performing  one  time  step  should  be  low  (0(1)). 
Interactions  should  be  local,  so  that  each  particle  directly  interacts  with  only  a  small 
neighborhood.  This  neighborhood  can  be  defined  based  on  the  topology  of  the  trian¬ 
gulation.  Also  the  convergence  should  be  fast.  A  minimal  amounts  of  time  steps 
should  be  taken  before  the  system  reaches  equilibrium. 

•  Robustness.  The  algorithm  should  always  converge. 

•  Quality.  The  resulting  mesh  should  locally  match  the  desired  edge  length.  For  now 
we  will  consider  this  length  equal  to  d,  which  is  defined  by  the  user.  In  a  more 
sophisticated  scheme,  the  implicit  geometric  feature  size  [25]  could  also  be  considered. 

The  chapter  is  structured  as  follows.  Section  3.1  introduces  the  general  dynamics  of  the 
interaction  model.  Section  3.2  defines  the  interaction  neighborhood.  Section  3.3  presents 
different  potential  functions.  Section  3.4  defines  a  model  for  controlling  the  particle  popu¬ 
lation.  Section  3.5  presents  some  improvements  regarding  speed.  Finally,  results  are  shown 
in  Section  3.6. 


3.1  Interaction  Model 

Whereas  in  Nature  and  in  some  previously  implemented  models  [27,  29],  the  motion  of  the 
particles  is  defined  by  a  second  order  differential  equation  {mx  =  F^{x,x)),  we  for  now 

consider  a  simpler,  first  order  model,  as  used  in  [34].  Although  first  order  models  are  more 
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likely  to  get  stuck  in  local  minima,  they  are  generally  numerically  more  stable.  The  motion 
of  the  particles  can  be  expressed  as 

(IT 

=  S7,E  +  C  (3.1) 

at 

where  x  is  the  2n-vector  containing  the  2-dimensional  positions  of  all  n  particles,  VxE  the 
gradient  of  the  potential  energy  field,  and  C  the  set  of  forces  imposed  by  domain  boundary 
conditions. 

The  behavior  of  the  system  highly  depends  on  the  choice  of  the  energy  potential  function 
E,  which  can  be  expressed  as  the  sum  of  the  potential  energies  between  each  pair  of  particles: 

E  =  Y.T.^''  (3.2) 

«  i 

where  E^^  is  the  potential  energy  of  particle  i  due  to  particle  j.  E^^  =  la),  where 

$  is  the  potential  function,  depends  only  on  the  distance  between  particles  i  and  j,  and 
the  desired  edge  length  a,  thus 

Vx^E^^  =  0  if  i  A:  and  j  A:  (3.3) 

and 

=  -I/*'  (t) 

where  —  a;*  is  defined  as  the  displacement  vector  between  particles  i  and  j,  a;*  being 

the  position  of  particle  k. 

Symmetry  is  also  a  desirable  feature,  thus 

(3.5) 

Finally,  interactions  are  local.  Particle  i  only  interacts  with  a  neighborhood  N\  thus 

=  0  if  j  ^  (3.6) 

Now  that  the  properties  of  the  energy  functions  are  better  defined,  it  is  possible  to  say 
more  about  its  gradient: 

Vx^E  = 

j  k 

= 

i  k 

i  k 

j 

and  since 

Vx^E  =  2-Y,\Ix^E^^  (3.7) 

j 
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(3.8) 


Also,  particle  i  only  interacts  with  particles  within  a  neighborhood  N^: 

V..E  =  2-Y,  V..E«  =  -2  i:  Jjf'  (^)  r’J 

jeN'  jeN'  ^  ^ 

To  simplify  the  above  expression,  we  define  A*-^  =  /a  as  the  “normalized”  distance 
between  particles  i  and  j: 


V,.T;  =  -2  ^  (3.9) 

jeN' 

Several  potential  functions  $  are  discussed  in  Section  3.3. 

3.1.1  Boundary  Conditions 

Boundary  condition  forces  are  necessary  to  make  sure  that  no  particle  moves  outside  the 
domain.  When  a  particle  i  is  on  a  constrained  vertex,  its  motion  should  be  zero,  thus 
C*  =  —Vx'E-  If  it  lies  on  a  domain  boundary,  motion  should  be  constrained  to  the  edge  it 
lies  on,  thus  C*  =  —  (n  •  V ^'E)  ■  n  where  n  is  a  unit  vector  normal  to  the  boundary.  For  all 
other  particles  C*  =  0. 

3.1.2  Numerical  Resolution 

To  solve  the  motion  equation  3.1,  numerical  integration  can  be  achieved  with  Euler’s  method: 

x{t  +  \t)  =  x{t)  +  ■  {V :,E  +  C)  (3.10) 

The  notion  of  time  can  be  eliminated  and  the  process  expressed  in  a  more  algorithmic 
fashion  as 


x^^x  +  5-{V:,E  +  C)  (3.11) 

Thus  for  each  particle  i  we  have 

x^  ^^x^ +5-{V^iE  +  &)  (3.12) 

3.1.3  Asynchronous  Updating 

Although  it  is  common  usage  to  update  the  positions  of  the  particles  synchronously,  we 
proceed  in  an  asynchronous  manner: 

loop 

randomly  pick  a  particle  i 
compute  V^iE  and  C* 
x^  ^  x^  +  5-ly^iE  +  C^) 
update  the  triangulation 

where  the  triangulation  update  step  is  performed  by  first  removing  site  i  and  then  reinserting 
it  at  its  new  position  (see  Section  2.5  on  moving  sites).  Most  of  the  time  the  topology  of  the 
mesh  remains  unchanged,  sometimes  the  motion  of  the  particle  can  result  in  the  swapping 
of  one  or  several  edges  in  the  triangulation. 

This  asynchronous  scheme  has  several  advantages  over  a  synchronous  one,  namely 
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•  it  introduces  some  noise  into  the  system  and  reduces  the  chances  of  reaching  a  local 
minimum 

•  it  opens  the  possibility  of  biasing  the  pick  and  accelerating  the  convergence  (more  on 
this  in  Section  3.5) 


3.1.4  Second  Order  Model 

Although  first  order  integration  is  nice  and  simple,  its  performance  can  be  poor,  as  the 
system  often  converges  to  a  not  so  good  local  minimum.  To  overcome  this  problem,  a 
second  order  model  is  introduced.  Its  formulation  is  rather  non-standard  but  it  is  easy  to 
implement,  and  doesn’t  seem  to  suffer  from  numerical  instabilities. 

u*  t—  13-v^  +  VxE  +  C 

a;*  t—  a;*  +  (5  •  u*  v  •  1 

The  velocity  vector  u*  can  be  seen  as  a  buffer  which  accumulates  the  forces,  and  which 
decays  exponentially.  The  differential  equation  associated  with  this  model  is  of  the  form 


X  +  k  •  X  =  VxE  +  C 


3.2  Interaction  Neighborhoods 

In  previous  work  [34,  27,  30,  29]  the  interaction  neighborhood  of  a  particle  i  was  defined 
based  on  geometrical  properties.  The  neighborhood  was  defined  by  a  sphere  with  a  radius 
equal  to  some  constant  times  the  desired  edge  length.  The  set  of  particles  77*  interacting 
with  a  particle  i  was  expressed  as 

Ka  =  {i  I  *  7^  h  <  cd} 

where  d*-^  is  the  distance  between  particles  i  and  j,  and  c  is  a  constant,  typically  between 
1.5  and  2.  For  fast  location  of  the  particles  within  77**,  k  —  d  trees  [29],  buckets  or  some 
other  additional  data  structures  are  used.  In  the  present  case,  however,  a  triangulation  is 
already  provided,  and  it  is  thus  easy  to  define  a  neighborhood  based  on  topology  instead 
of  geometry.  The  distance  between  two  particles  i  and  j  is  then  defined  by  the  number  of 
edges  in  the  shortest  path  between  i  and  j  in  the  triangulation.  The  set  of  particles  77* 
interacting  with  a  particle  i  becomes 

—  {j  \  ^  h  3path  between  i  and  j  of  length  <  k} 

In  practice  small  neighborhoods  are  desirable,  so  only  path  lengths  of  A:  =  1,2  are 
considered.  Using  the  quadedge  terminology  from  Chapter  2,  these  neighborhoods  are  (also 


see  Figure  3.1): 

nI 

=  0 

3e  :  e.Org  =  i,  e.Dest  =  j} 

(3.14) 

Ni 

=  0 

\il^  3,  3e,  /  :  e.Org  =  i,  e.Dest  =  f.Org,  f.Dest  =  j} 

(3.15) 

The  size  of  the  neighborhoods  are  typically  6  for  77i,  and  18  for  772.  Since  they  are  both 
0(1)  in  size,  an  update  operation  as  described  in  Section  3.1.3  requires  only  0(1)  time, 
whereas  it  would  take  0(n)  if  all  particles  interacted  with  all  others. 
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Ni 


N2 


Figure  3.1:  Two  different  sizes  of  neighborhoods 


It  is  of  course  also  possible  to  define  a  neighborhood  based  on  geometrical  and  topological 
properties. 


3.3  Potential  Functions 

This  section  presents  several  potential  functions  that  have  been  used  in  the  past,  but  also 
introduces  new  ones.  All  potential  functions  should  fulfill  the  following  desirable  properties: 

•  continuity.  The  function  and  its  first  derivative  should  be  continuous  functions 

•  no  singularities.  The  function  should  have  finite  values  for  any  finite  inter-particle 
distance. 

3.3.1  Requirements  at  Equilibrium 

The  optimal  particle  arrangement  is  a  regular  array  of  equilateral  triangles  with  edge  length 
of  (7.  The  forces  between  particles  should  be  defined  so  that  this  pattern  is  an  equilibrium. 
Moreover  the  equilibrium  should  be  stable,  that  is,  if  the  positions  of  the  particles  are  slightly 
altered,  they  should  move  back  to  the  optimum. 

These  considerations  should  also  hold  for  particles  near  or  on  the  boundaries  of  the 
domain.  If  particles  on  the  boundaries  are  constrained  to  stay  on  it,  no  problem  arises.  But 
particles  which  are  one  edge  length  away  from  the  boundary  are  not  at  equilibrium  if  the 
neighborhood  is  N2  or  larger:  if  the  particles  repel  each  other  at  distances  larger  than  a, 
they  tend  to  drift  towards  the  boundary,  and  in  the  reverse  direction  if  the  interaction  is 
attractive.  Thus,  to  guarantee  optimality  near  the  boundaries  the  interaction  neighborhood 
should  be  . 

On  the  other  hand  N2  neighborhoods  are  more  powerful  at  organizing  particles  away 
from  the  boundaries.  A  tradeoff  must  thus  be  found  between  quality  near  the  boundaries 
and  elsewhere.  As  the  effects  depend  a  lot  on  the  choice  of  the  potential  energy  function, 
we  empirically  search  for  the  best  solution  for  each  function. 


31 


3.3.2  Gaussian  Potential 

Witkin  and  Heckbert  [34]  used  a  Gaussian  potential  function  to  spread  particles  on  implicit 
surfaces  (Figure  3.2): 

=  exp  (-2(A*-»')2) 


Potential  $(A*'^)  Derivative  $'(A*'^) 


Figure  3.2:  Gaussian  potential  function 

This  model  works  well  with  N2  neighborhoods  but  fails  to  arrange  particles  into  the 
desired  pattern  with  N^.  Originally  Witkin  and  Heckbert  used  neighborhoods. 

Instead  of  having  the  potential  be  a  Gaussian  function,  it  is  possible  to  have  it  be  the 
derivative  of  a  Gaussian  (Figure  3.3): 

$(A*-»')  =  A*-»'exp  (3.16) 

of  which  the  derivative  is 

$'(A*^')  =  (1  -  (A*^)2)exp  (3.17) 


Potential  $(A*'^)  Derivative  $'(A*'^) 


Figure  3.3:  Gaussian  derivative  potential  function 
This  new  scheme  generates  regular  patterns  for  Ni  and  N2  neighborhoods. 
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3.3.3  Lennard-Jones  Potential 

One  of  the  most  widely  used  interaction  model  [29,  27]  is  based  on  the  Lennard-Jones 
potential,  or  van  der  Waals  force  which  is  inspired  from  molecular  chemistry.  The  potential 
is  defined  as  (Figure  3.4) 

-  6(A*-»y2  (3.18) 


0  0.2  0.4  0.6  0.8  1  1 .2  1 .4  1.6  1 .8  2 


Potential  $(A*'^) 


Derivative  $'(A*'^) 


Figure  3.4:  Lennard-Jones  potential  function 

The  problem  with  such  a  model  is  the  singularity  at  A*-^  =  0  ($'(0+)  oo).  Numerical 
integration  can  thus  be  disastrous. 

Shimada  [27]  has  defined  a  potential  function  with  a  similar  shape  but  which  overcomes 
the  singularity  problem  (Figure  3.5): 

$'(A*^')  =  ^  (A«y  -  ^  (A«y  +  5  (3.19) 


Potential  $(A*'^)  Derivative  $'(A*'^) 


Figure  3.5:  Shimada’s  potential  function 

Shimada  used  but  this  scheme  also  works  with  Ni  and  N^. 

Other  approximations  of  the  Lennard-Jones  potential  function  are  possible,  such  as  (Fi¬ 
gure  3.6) 

$'(A*-^')  =  (1  -  (A*-»y)  exp(-d(A*-»y)  (3.20) 
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1.2 


Derivative  $'(A*'^) 


Figure  3.6:  New  Lennard-Jones  approximation 


3.3.4  Laplacian  Smoothing 

Laplacian  smoothing  [8]  is  a  widely  used  technique  for  improving  the  shape  of  triangles  in 
a  mesh.  In  particle  terminology,  Laplacian  smoothing  consists  of  moving  a  particle  i  to  the 
centroid  of  its  neighbors.  Its  displacement  can  be  expressed  as 


Ax^  = 


jeNi 


jdNl 


As  the  displacement  is  proportional  to  the  derivative  of  the  potential,  the  latter  can  be 
defined  as  (Figure  3.7) 


(3.21) 


Potential  $(A*'^)  Derivative  $'(A*'^) 


Figure  3.7:  Laplacian  potential  function 


3.3.5  Error  Potentials 

Another  way  of  putting  the  problem  is  to  define  a  quality  measure  of  the  mesh  and  consider 
the  particle  interaction  as  an  optimization  process  by  gradient  descent. 
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We  here  define  two  quality  measures.  The  first  one  is  based  on  the  difference  between 
the  actual  and  the  optimal  distances  between  particles  (Figure  3.8): 

$(A«)  = -i(Ay  -  1)2  (3.22) 

of  which  the  derivative  is 

$'(A«)  =  1  -  A*-^'  (3.23) 


Potential  $(A*'^)  Derivative  $'(A*'^) 


Figure  3.8:  Error  potential  function 

The  second  quality  measure  is  based  on  the  ratio  of  the  distances  (Figure  3.9): 

$(Ab)  =  -ilog2A«  (3.24) 

Although  this  definition  can  seem  obscure  at  first,  one  can  easily  verify  that  error  is  the 
same  for  rfP  =  P  ■  a  and  ^  •  d,  that  is,  edges  that  are  too  long  or  too  short  by  a  factor 

P  are  equally  penalized.  The  derivative  of  the  potential  is: 

$'(A*^')  =  (3.25) 
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Potential  $(A*'^) 


Derivative  $'(A*'^) 


Figure  3.9:  Error  potential  function 
The  natural  interaction  neighborhood  of  these  potentials  is  Ni . 
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3.3.6  Which  is  best? 

Empirical  testing  has  led  to  the  following  conclusions: 

•  models  that  are  attractive  and  repulsive  perform  better  than  the  ones  that  are  attra¬ 
ctive  only  or  repulsive  only 

•  the  N2  neighborhood  does  a  better  job  at  arranging  particles  to  create  equilateral 
triangle  away  from  the  boundaries. . . 

•  . . .  but  performs  rather  poorly  near  the  boundaries 

In  practice  a  Lennard-Jones  like  scheme,  such  as  the  proposed  approximations,  seems 
best. 


3.4  Adaptive  Population  Control 

Generally  it  is  not  known  a  priori  how  many  particles  are  needed  for  the  edges  to  conform 
the  desired  size.  The  number  could  be  determined  from  the  area  of  the  domain,  but  there 
are  several  disadvantages  to  this  strategy: 

•  when  generalizing  to  surfaces  in  3-D  or  to  a  non-uniform  desired  edge  length,  it  is 
difficult  to  correctly  estimate  the  area. 

•  even  if  one  knew  how  many  particles  are  needed  to  optimally  fit  the  domain,  the 
question  would  remain  of  where  to  put  them.  If  the  initial  configuration  is  poor,  it 
can  take  a  long  time  for  the  system  to  reach  a  reasonably  good  configuration. 

An  alternative  strategy  is  an  adaptive  scheme  where  new  particles  are  created  where 
their  density  is  too  low,  and  others  are  annihilated  where  too  high.  The  problem  thus 
becomes  to  evaluate  particle  density.  In  a  first  step  we  will  consider  the  1-dimensional  case, 
and  then  generalize  to  two  dimensions. 


3.4.1  1-D  Algorithm 

Let  5  be  a  set  of  edges,  and  L{S)  the  sum  of  the  normalized  lengths^  of  each  segment  in  S. 
Ideally  L{S)  should  be  equal  to  L{S)  =  |5|. 

The  particle  density  of  S  is  defined  as  S{S)  =  L{S)/L{S).  The  density  error  measure 
associated  with  S  is  defined  to  be 


f  S(S)  if  S(S)  >  1 
I  1/S{S)  otherwise 


(3.26) 


Creation  Rule 

Let  S  be  composed  of  one  edge  e,  and  S'  the  set  of  the  two  halves  of  e  after  a  particle  has 
been  inserted  at  its  middle  point.  A  particle  is  inserted  if  s{S')  <  s{S).  We  know  that 
L{S)  =  L{S')  =  ||e||,  |5|  =  1,  and  |5'|  =  2.  Thus  a  new  particle  should  be  inserted  if 
||e||  >  v^,  where  ||e||  is  the  normalized  length  of  edge  e. 

^the  normalized  length  is  defined  to  be  the  actual  length  divided  by  a 
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Annihilation  Rnle 


Let  S  be  the  two  edges  ei  and  62  on  each  side  of  a  particle  p.  Let  S'  be  a  set  containing  the 
concatenation  e  of  ei  and  62  after  removal  of  p.  Particle  p  is  annihilated  if  s{S')  <  s{S). 
We  know  that  L{S)  =  L{S')  =  ||ei||  +  ||e2||,  |5|  =  2,  |5'|  =  1.  Thus  particle  p  should  be 
annihilated  if  ||ei||  +  ||e2||  <  v^,  where  ||ej||  is  the  normalized  length  of  edge  e,. 


3.4.2  2-D  Algorithm 


Let  T  be  a  set  of  triangles,  and  A{T)  the  sum  of  the  normalized  areas^  of  each  triangle  in 
T.  Ideally  A{T)  should  be  equal  to  A{T)  =  |T|A,  where  A  =  -s/S/d  is  the  area  of  a  unit 
equilateral  triangle. 

The  particle  density  of  T  is  defined  as  S{T)  =  A{T)/A{T).  The  density  error  measure 
associated  with  T  is  defined  to  be 


f  S{T)  if  S{T)  >  1 
I  1/S{T)  otherwise 


(3.27) 


Creation  Rule 

Let  T  be  the  set  of  the  triangles  on  each  side  of  an  edge  e,  and  T'  the  set  of  triangles  incident 
to  the  two  halves  of  e  after  a  particle  has  been  inserted  at  its  middlepoint.  A  particle  is 
inserted  if  e(T')  <  e(T).  We  know  that  A{T)  =  A{T'),  |T|  =  2,  and  |T'|  =  4.  Thus  a  new 
particle  should  be  inserted  on  e  if  A{T)  >  a/3/2. 

A  more  sophisticated  rule  takes  into  consideration  the  local  optimization  of  the  topology 
(Delaunay  triangulation)  after  insertion.  The  average  degree  of  a  node  inside  a  mesh  is  6, 
and  |T'|  should  thus  be  6.  Since  it  would  be  costly  to  compute  the  areas  of  the  triangles 
after  the  local  optimization,  we  consider  two  virtual  triangles  which  are  affected  by  the  edge 
swapping.  Let  U  be  the  set  T  augmented  by  two  virtual  triangles  of  ideal  size,  and  U'  the 
set  U  augmented  by  the  two  same  virtual  triangles.  We  know  A{U)  =  A{U')  =  A{T)  +  2A, 
|17|  =  4,  and  |17'|  =  6.  Thus  a  new  particle  should  be  inserted  on  e  if  A{T)  >  3/v^  —  ■s/3/2. 

Annihilation  Rule 

Let  T  be  the  set  of  triangles  incident  to  a  particle  p 
from  the  retriangulation  of  the  region  covered  by  T 
e(T')  <  e(T).  We  know  that  A(T')  =  A{T)  and  |T'| 
annihilated  if  A{T)  <  a/3|T|(|T|  -  2)/4. 

3.4.3  Combination  of  1-D  and  2-D  rules 

The  rules  that  are  actually  used  in  our  algorithm  are  a  combination  of  1-dimensional  and 
2-dimensional  rules.  Let  i  be  the  particle  that  is  picked  at  the  given  step. 

First,  the  annihilation  rules  are  considered.  If  the  particle  i  lies  on  a  boundary,  then  the 
1-D  rule  in  applied.  Otherwise  the  2-D  rule  is  applied. 

If  the  particle  is  not  annihilated,  the  creation  rules  are  considered  for  all  the  edges 
that  are  connected  to  particle  i.  If  several  edges  are  candidates  for  insertion,  the  one  that 
maximizes  the  improvement  in  particle  density  is  selected  for  splitting. 

^the  normalized  area  is  defined  to  be  the  actual  area  divided  by 


.  Let  T'  be  set  of  triangles  resulting 
.  Clearly  p  should  be  annihilated  if 
=  |T|  —  2.  Thus  particle  p  should  be 
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3.4.4  Initial  Population 

The  initial  population  corresponds  to  the  sites  of  the  PSLG  defining  the  domain.  These 
particles  are  never  moved  or  annihilated  so  that  the  triangulation  always  remains  confor¬ 
ming. 


3.5  Speeding  up  the  Process 

Since  the  particles  are  moved  one  by  one,  it  is  possible  to  define  a  heuristic  so  that  the  ones 
that  “need”  to  be  moved  first  are  indeed  moved  first.  The  introduction  of  an  alive  flag  for 
each  particle  can  lead  to  a  simple  heuristic: 

•  all  the  particles  are  initially  alive 

•  when  an  alive  particle  p  is  picked  and  updated,  if  its  normalized  speed  |u^|/(T  is  larger 
than  a  threshold,  all  of  its  dead  neighbors  become  alive.  If  the  speed  is  too  low,  p  is 
marked  dead. 

•  when  a  particle  p  is  annihilated,  all  of  its  dead  neighbors  become  alive. 

•  when  a  particle  p  is  created,  p  and  all  of  its  dead  neighbors  become  alive. 

•  a  dead  particle  is  never  picked.  Consequently,  particles  don’t  move  while  they  are 
dead. 

More  complex  schemes  have  been  tested  such  as  assigning  picking  probabilities  to  each 
particle  proportional  to  their  speed  or  their  total  energy  (kinetic  and  potential) .  Preliminary 
results  have  shown  no  improvement  over  the  simple  scheme.  Moreover  the  computational 
cost  of  such  a  method  is  high  (time  cost  for  one  step  would  be  O(logn)  instead  of  0(1)). 

3.5.1  Ending  the  simulation 

The  heuristic  defined  in  the  previous  section  allows  a  very  simple  test  to  end  the  simulation: 
end  when  no  more  particles  are  alive. 


3.6  Results 

In  this  section  we  briefly  present  results  for  a  simple  mesh.  The  domain  that  is  meshed  is 
a  unit  square,  and  the  desired  edge  length  is  set  to  a  =  0.02.  The  mesh  we  have  obtained 
with  an  approximated  Lennard-Jones  potential  (Equation  3.20)  is  depicted  in  Figure  3.10. 

The  angle  and  normalized  edge  length  histograms  show  that  the  elements  are  close  to 
optimality  (i.e.  60  degree  angles,  and  unit  normalized  edge  length). 
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Normalized  edge  length  histogram 


Figure  3.10:  A  simple  mesh 
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Chapter  4 

Interlude:  Delaunay 
Triangulations  with  Java 


Many  people  are  unfamiliar  with  computational  geometry  and  Delaunay  triangulations.  The 
World  Wide  Web  appears  to  be  a  good  tool  to  familiarize  those  people  with  these  concepts. 
And  as  learning  is  often  easier  when  playing  games,  we  have  designed  an  interactive  Delaunay 
triangulator  using  the  Java  programming  language  developed  by  Sun  [12].  It  allows  all  people 
using  a  Java  enabled  browser  such  as  Netscape  Navigator'^“  2.0  to  build  a  triangulation  by 
interactively  inserting  and  removing  sites.  The  triangulator  can  be  found  at: 

http : //www . cs . emu . edu/~bossen/ triangulator . html 

Complete  user  instructions  are  given  on  the  web  page. 

The  algorithm  that  has  been  implemented  is  the  incremental  one.  When  the  user  adds 
a  new  point,  it  can  simply  be  inserted  into  the  triangulation.  When  a  point  is  removed  the 
whole  triangulation  is  rebuilt^.  The  quadedge  data  structure  has  been  used  to  represent  the 
triangulation.  This  choice  is  justified  by  the  fact  that  the  Voronoi  diagram  is  also  computed 
and  displayed. 


^Although  not  very  effleient,  this  solution  was  easier  to  implement 
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Chapter  5 

Anisotropic  Meshes 


In  this  chapter  we  generalize  the  particle-based  mesh  generation  method  described  in  Cha¬ 
pter  3  to  nonuniform,  isotropic  and  anisotropic,  meshes.  The  feature  size  is  no  longer  a 
constant  d,  but  is  position-dependent  and  direction-dependent. 

In  strongly  directional  phenomena  such  as  shocks  and  limit  layers  in  fluid  flows,  the 
elements  should  not  only  adapt  in  size  but  also  in  shape.  They  should  be  stretched  in  the 
direction  of  the  flow.  The  desired  edge  length  depends  on  its  orientation,  hence  the  name 
anisotropy.  A  metric  tensor  is  introduced  to  quantify  the  stretching  of  the  triangles  at  every 
point  in  the  domain.  It  reflects  the  curvature  of  the  function  that  is  approximated. 

Relatively  little  work  has  been  done  on  anisotropic  meshes.  D’Azevedo  [5]  has  proposed 
the  use  of  coordinate  transformations  to  generate  optimal  triangulations.  Although  very 
elegant,  his  method  is  limited  to  metric  tensors  that  represent  a  fiat  space  (the  Riemann- 
Christoffel  tensor  [28],  that  is  the  curvature  of  the  space  represented  by  the  metric,  should 
be  zero).  Another  weakness  of  his  method  is  poor  meshing  at  boundaries. 

Peraire  et  al.  [23],  and  later  Mavriplis  [22]  introduced  stretching  vectors  to  quantify  the 
anisotropy,  which  really  are  eigenvectors  of  the  metric  tensor,  but  never  formally  introduced 
the  latter.  Greedy  insertion  algorithms  have  been  proposed  by  researchers  at  INRIA  [3,  31]. 
We  believe  that  better  node  positioning  can  be  achieved. 

This  chapter  is  structured  as  follows.  Section  5.1  introduces  some  concepts  of  Rieman- 
nian  geometry.  Section  5.2  exposes  the  changes  made  to  the  triangulation  procedure  due 
to  anisotropy.  Section  5.3  presents  an  anisotropic  version  of  the  particle  interaction  mo¬ 
del  described  earlier.  Section  5.5  summarizes  the  meshing  algorithm.  Section  5.4  presents 
background  meshes  which  are  convenient  for  representing  the  metric.  Section  5.6  shows  how 
to  derive  the  metric  M  for  a  given  function  /  that  the  mesh  should  adapt  to.  Section  5.7 
shows  how  the  mesh  can  incrementally  be  adapted  to  a  function  when  no  a  priori  knowledge 
of  the  solution  is  available. 


5.1  Riemannian  Geometry 

In  Euclidean  geometry,  the  definition  of  distances  is  isotropic  and  very  simple.  Riemannian 
geometry  describes  “warped”  or  curved  spaces.  In  Riemannian  geometry,  the  anisotropy  of 
distance  is  defined  by  a  metric  tensor  M.  This  tensor  quantifies  the  desired  stretching  of 
the  triangles  in  the  mesh.  M  is  a  symmetric  2  by  2  matrix,  and  both  of  its  eigenvalues  Ai 
and  A2  are  positive: 
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M  = 


(5.1) 


/  cos  0  sin  0  \  /  Ai  0  ^  ^  —  sin  0  A 

—  sin  0  cos  0  J  0  A2  J  y  sin  0  cos  6  J 

If  both  the  eigenvalues  are  equal  to  1,  then  M  =  I  and  the  space  is  Euclidean.  If  M  is 
considered  locally  constant,  the  unit  ball  is  the  ellipse  px'^Mx  =  1.  The  first  axis  of  the 
ellipse  has  length  n  =  I/1/A1,  and  makes  an  angle  9  with  any  horizontal  line.  The  second 
axis  has  length  r2  =  (Figure  5.1). 


Figure  5.1:  Elliptic  unit  ball 

In  Riemannian  geometry,  the  basic  geometrical  operators  are  redefined.  The  dot  product 
becomes 

a^M  b  (5.2) 

and  the  cross  product 

V det  M {a  x  b)  (5.3) 

where  a  x  b  =  aib^  —  a^bi  is  the  2-dimensional  cross  product. 

5.1.1  Computing  Distances 

In  Riemannian  geometry,  the  length  of  a  parametric  curve  T(t)  between  points  i  and  j, 
where  t  €  [0, 1],  r(0)  =  x\  and  r(l)  =  x^  is  defined  as 

m  =  f  ^f(t)^M(r(t))f(t)  dt  (5.4) 

where  h  =  ^.  The  distance  between  two  points  is  the  length  of  the  shortest  parametric 
curve  r.  Such  a  curve  is  called  a  geodesic.  Computing  the  geodesic  for  an  arbitrary  metric 
is  not  trivial.  It  is  simpler  to  approximate  the  distance  by  computing  the  length  of  a  simple, 
non-guaranteed  shortest,  path.  Integrating  along  the  straight  line  T{t)  =  x^  +  t  ■  {x^  —  a;*) 
seems  reasonable.  If  M  varies  linearly  between  i  and  j  then  the  distance  between  the  two 
points  can  be  defined  as 

Y^rd^(M*  +  t{M^  —  M*))r*-1  dt 

^{d^y  +  t{{d^y-{d^y)  dt  (5.5) 
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where  is  the  distance  as  “seen”  from  point  k.  With  appropriate  variable 

changes,  the  integral  solves  to 


2(d^)^+d^d^  +  (d^r 
3  df  +  4^ 


(5.6) 


Assuming  that  particles  i  and  j  are  close  to  each  other,  and  that  the  metric  is  slowly 
changing,  cruder  approximations  of  the  distance  are  also  possible,  namely  by  simply  avera¬ 
ging  the  distances  as  “seen”  from  i  and  j 


or  by  averaging  the  metrics  as 


(5.7) 


- ^-T  •  • 

K  V  ^  with  ^ -  (5.8) 

From  the  three  proposed  approximations  the  last  is  the  fastest  to  compute  since  it  requires 
only  one  square  root  operation.  Also  finding  the  squared  distance  doesn’t  require  any  square 
root  operation  at  all. 


5.1.2  Computing  Areas 

The  area  of  a  domain  is  defined  as 


A(ff) 


\/ det M (x)  dxidx2 


(5.9) 


JJn 

Again  it  would  be  possible  to  solve  the  integral  for  a  triangle  Aabc,  knowing  the  metric 
at  each  vertex,  and  linearly  interpolating  it  in  between. 

We  will  not  bother  to  do  so,  and  simply  give  a  cruder  approximation: 


A{Aabc) 


+M’>  + 


{b  —  a)  X  {c  —  a) 


(5.10) 


5.2  Anisotropic  Triangulation 

Although  it  is  possible  to  compute  the  Delaunay  triangulation  of  an  anisotropically  distri¬ 
buted  point  set,  it  is  reasonable  to  consider  anisotropicity  in  the  triangulation  process.  It 
has  been  shown  by  D’Azevedo  [5]  that  it  is  better  to  Delaunay  triangulate  in  a  transformed 
space  to  minimize  errors  for  function  interpolation. 

Only  simple  cases  (for  example  when  the  metric  tensor  is  constant)  allow  a  global  map 
of  the  domain  into  a  2-dimensional  Euclidean  space  by  coordinate  transformation.  In  the 
general  case,  spaces  of  higher  dimensions  are  needed.  To  solve  this  problem,  a  simple 
approximation  is  used.  It  is  considered  that  the  metric  is  locally  constant  when  performing 
the  circum circle  test.  This  local  constant  metric  is  computed  by  averaging  the  metric  at  the 
four  points  considered  for  the  test: 

1  4 

i=l 
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In  Euclidean  geometry  the  InCircle  test  is  defined  by  equation  2.9.  Its  Riemannian 
equivalent  is: 

a/ det  M (a  x  h){(F M  d)  +  {a^M  b)  V det  M(c  x  d)  >  0  (5-12) 

and  since  Vdet  M  is  positive: 

{a  X  b){c^ M  d)  +  {a^ M  b){c  X  d)  >  0  (5.13) 


5.3  Anisotropic  Energy  Potential 

All  the  energy  potential  functions  defined  in  Chapter  3  are  in  the  form 


and  the  gradient  of  which  is 


\ii 


(5.14) 


In  an  anisotropic  context  the  normalized  distance  is  replaced  by  d^  as  defined  in 
Section  5.1.1.  Finding  an  equivalent  for  is  more  difficult.  In  it  could  be  replaced  by 
l/V det  M* ,  but  that  would  break  the  symmetry  between  and  E^^ .  To  solve  this,  we 
can  average  1  /  V det  M*  and  1  /  V det  to  get  .  As  this  expression  is  rather 

complex,  in  practice  we  will  approximate  it  by  /t  =  l/V det  Mb' ,  where  is  the  average 
between  M*  and  ME 

The  potential  energy  function  becomes 


E^^  =  K  ■  $(d*-»') 


of  which  the  gradient  is 


V^iE^^  =  K  ■  ^'{d^^)  ■  V^id^^  +  $(d*-»')  •  V^iK 

where  is  quite  a  complex  expression,  which  depends  on  the  derivative  of  the  position 

of  particle  i,  but  also  on  the  derivatives  of  M.  For  the  sake  of  simplicity,  we  will  consider  that 
^  as  defined  is  equation  5.8,  and  that  is  locally  constant.  Following 
this  simplification  we  have 


Mb'rb 

db' 


Furthermore  is  zero  since  is  considered  locally  constant.  The  gradient  can 

thus  be  rewritten  as 


VriE^^ 


$'(d*J) 

db' 


(5.15) 


It  is  also  possible  to  simplify  the  expression  of  the  gradient  in  a  more  radical  way,  by 
substituting  d*-^  for  A*-^  in  the  gradient  expression  5.14: 


db' 


(5.16) 
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Equations  5.15  and  5.16  are  quite  similar,  but  whereas  in  the  latter  the  gradient  of 
is  aligned  with  it  is  usually  not  the  case  in  the  former,  where  is  multiplied  by  the 
matrix  This  multiplication  has  for  consequence  a  precession  movement  that  will 

tend  to  align  with  one  of  the  eigenvectors  of  M’-E  When  is  parallel  to  one  of  the 
eigenvectors  of  then  the  gradient  of  is  aligned  with  . 

In  practice,  the  simpler  scheme  (Equation  5.16)  is  used. 


5.4  Background  Mesh 


Although  it  is  possible  to  explicitly  define  the  metric  tensor  as  a  set  of  functions,  it  is  often 
preferable  to  define  the  metric  at  given  points,  and  then  interpolate  in  between.  A  backgro¬ 
und  mesh  is  used  to  do  so.  At  each  node  of  the  background  mesh  a  metric  matrix  is  defined. 
The  metric  at  the  other  points  of  the  domain  can  be  computed  by  linear  interpolation.  If 
M“,  and  are  the  metric  matrices  at  the  vertices  of  a  triangle  l^abc  in  which  a  point 
p  lies,  then  the  metric  matrix  at  p  is  given  by: 


y^ap  .  ^  ^bp  .  j^b  ^  ^cp  . 

MP  =  - r - 

yjap  q.  yjbp  q.  yjcp 


where 


s 

II 

(xP 

-x’>) 

X 

{x^ 

-x'>: 

s 

II 

{xP 

-x^) 

X 

(a;“ 

-x^ 

W^P  = 

(xP 

-x^) 

1  X 

(x'> 

-x^ 

(5.17) 


The  triangulation  of  the  set  of  samples  defining  the  metric  can  either  be  given  by  the 
user  or  computed.  In  the  latter  case,  Delaunay  triangulation  is  used. 

To  compute  the  metric  at  a  given  point,  the  triangle  it  lies  in  first  has  to  be  determined. 
To  do  so,  it  is  possible  to  use  data  structures  such  as  the  history  tree  proposed  in  [14]. 
Although  the  cost  of  one  location  is  only  0(log  to)  ,  where  to  is  the  number  of  samples  defining 
the  metric,  this  method  is  suboptimal.  Since  the  metric  is  only  evaluated  at  positions  where 
particles  are  located,  and  that  particles  rarely  move  across  more  than  one  triangle  in  the 
background  mesh  during  one  time  step,  it  is  better  to  cache  the  triangle  a  particle  lies 
in.  When  a  particle  is  moved,  the  eventually  new  triangle  it  moves  into,  can  quickly  be 
determined  by  a  simple  walking  method.  The  metric  evaluation  cost  can  thus  be  reduced 
to  0(1). 

A  potential  drawback  of  only  piecewise  linearly  defining  the  metric  is  that  the  feature 
size  function  only  has  O®  continuity.  In  the  present  scheme,  however,  this  is  not  a  problem. 


5.5  Method  Summary 

Section  3.1.3  introduced  the  main  loop  of  the  mesh  generator.  We  will  now  extend  it  to  the 
general,  anisotropic  case: 

build  the  constrained  Delaunay  triangulation  of  the  domain 

loop 

randomly  pick  a  particle  i  that  is  marked  alive 
compute  Vx'E  and  C* 
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^  +  5-{V^iE  +  C^) 

compute  the  areas  of  triangles  incident  to  i 
create/annihilate  a  particle  if  areas  are  too  high/low 
update  the  triangulation 

where  equations  5.15  and  3.20  are  used  to  compute  and  equation  5.10  to  compute 

triangle  areas.  For  distances,  either  equation  5.8  or  5.6  can  be  used. 

This  algorithm  has  been  implemented  in  C++  on  Silicon  Graphics  workstations.  The 
source  code  spans  over  roughly  5000  lines.  The  OpenGL  graphical  library  was  used  to 
display  the  mesh  evolution  in  real  time.  All  the  meshes  inserted  in  this  document  have  been 
generated  by  the  program,  which  can  directly  output  encapsulated  Postscript  files. 


5.6  Function  Interpolation 


Given  a  function  /,  the  metric  that  generates  a  mesh  for  piecewise  linear  interpolation  such 
that  the  RMS  error  is  minimum,  is  given  by  the  “curvature”  of  the  function.  This  curvature 
is  equivalent  to  the  Hessian  matrix  of  the  function  defined  as 


(  d'^f/dxj  d‘^f/dxidx2  \ 
V  d‘^fldx2dxi  d'^f/dxl  ) 


Since  E[^  is  symmetric,  it  can  be  represented  as 


Rf  = 


Ai  0 

0  A2 


R 


(5.18) 


(5.19) 


where  i?  is  a  rotation  matrix.  Given  this  decomposition,  it  is  trivial  to  define  a  matrix 
such  that  both  its  eigenvalues  are  positive^ : 

'o'  |a",|)«  <“<» 

In  practice  it  is  necessary  to  limit  the  values  of  the  eigenvalues  to  a  given  range.  Other¬ 
wise  the  mesh  elements  could  get  arbitrary  small  or  large.  It  has  also  been  proposed  [3]  to 
align  one  of  the  eigenvectors  of  the  metric  matrix  with  the  boundary,  when  close  to  one. 


5.6.1  Example:  a  Gaussian  function 

The  function  we  approximate  is 


f{x,y)  =  exp 


of  which  the  Hessian  is 


Rf 


-  1  xy  \ 

xy  -  1  y  •' 


The  eigenvalues  of  R-^  are  Ai  : 
eigenvectors  are  ei  =  (  —yfr  xjr 


f  and  A2 
and  62  = 


=  (a;2  +  —  1)  •  /.  The  corresponding 

'j'  _ _ 

(  xjr  y/r  )  ,  where  r  =  ^/x‘^  +  y"^. 


^The  eigenvalues  need  to  be  positive  to  ensure  that  the  squared  distance  between  any  two  points  is  always 
positive. 
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Following  the  rules  previously  described,  we  obtain  the  metric  tensor 


f  y‘^  +  cx^  (c  -  l)xy  \  ^  |/| 

(c  -  l)xy  x"^  +  cy"^  )  x"^  +y‘^ 


where  c  =  \x‘^  +  y"^  —  1\. 

Figure  5.2  shows  the  mesh  obtained  with  the  over  the  domain  fl  =  [0,3]  x  [0,3]. 


Unadapted:  572  particles 


Adapted:  409  particles 


Figure  5.2:  Meshes  for  Gaussian  function  approximation 

Experiments  have  shown  that,  for  an  identical  number  of  particles,  the  RMS  error  of  the 
piece-wise  linear  approximated  function  is  about  8  times  smaller  for  adapted  meshes  than 
for  unadapted  ones. 


5.7  Incremental  Mesh  Adaptation 

To  generate  a  good  mesh  for  finite  element  analysis,  one  has  to  have  a  good  knowledge 
of  what  the  solution  is.  But  what  if  one  doesn’t?  A  solution  is  to  generate  a  metric  by 
successive  refinements.  The  incremental  adaptation  loop  is 

1.  start  with  a  uniform  feature  size  function 

2.  build/update  the  mesh 

3.  run  the  finite  element  solver  with  the  current  mesh 

4.  estimate  the  error.  If  the  error  is  smaller  than  some  amount,  then  exit,  else  adapt  the 
metric  according  to  the  current  solution,  and  go  back  to  step  2 

While  the  initial  mesh  is  uniform  as  in  Figure  5.3,  after  several  iterations,  the  metric 
becomes  better  adapted  to  the  problem  and  the  quality  of  the  mesh  is  improved,  similarly 
to  Figure  5.4^. 

^This  figure  has  been  obtained  with  a  hand-defined  background  mesh.  No  solver  has  been  used. 
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Figure  5.3:  Without  background  mesh:  uniform  mesh  with  2865  particles 


Figure  5.4:  With  background  mesh:  anisotropic  mesh  with  229  particles 


5.8  Running  Time 

To  empirically  determine  the  asymptotic  running  time  of  the  algorithm,  we  have  generated 
meshes  for  Gaussian  function  approximation  with  different  numbers  of  particles. 

Figure  5.6a  shows  the  CPU  time^  required  to  generate  a  mesh  containing  n  particles. 
The  running  time  increases  slightly  faster  than  linearly  with  the  number  of  particles.  To  find 
out  why  this  is,  we  have  looked  at  the  growth  of  the  number  of  steps  versus  the  number  of 
particles.  Figure  5.6b  shows  the  number  of  steps  grows  less  than  linearly  with  the  number  of 
particles.  On  the  other  hand,  the  number  of  steps  executed  within  a  second  decreases  as  the 
number  of  particles  grows.  This  is  in  contradiction  with  the  theory,  which  predicts  that  the 
time  complexity  of  one  step  is  0(1).  Our  guess  is  that  speed  performance  decreases  because 
the  number  of  cache  misses  grows.  Using  profiling,  we  have  found  that  the  subroutine  which 
gets  more  and  more  time  consuming  proportionally  to  the  number  of  steps,  is  a  quite  simple 
subroutine  which  does  actually  only  scan  through  all  the  edges  adjacent  to  a  newly  picked 

®On  a  250Mhz  Mips  R4400  chip 


50 


Figure  5.5:  Background  mesh 


Figure  5.6:  Running  time  results 

particle  i.  The  reason  for  apparent  sup-linear  time  complexity  is  thus  cache  misses. 
Our  conclusion  is  that  the  asymptotic  time  complexity  of  the  algorithm  is  0{n). 
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Chapter  6 

Conclusion 


We  have  presented  a  new  method  for  generating  triangular  meshes  in  a  2-dimensional  domain 
bounded  by  straight  lines.  The  following  considerations  have  led  to  the  formulation  we  have 
given: 

•  the  particle  approach  seemed  promising  because  it  allowed  good  node  placement. 

•  whereas  in  previous  methods  the  particle  interaction  neighborhood  was  defined  by 
geometrical  distances,  we  have  opted  for  topological  distances.  This  allows  for  faster 
neighborhood  computation  at  zero  data  structure  space  cost. 

•  although  it  is  possible  to  define  particle  interaction  as  purely  attractive  or  repulsive, 
it  turned  out  that  mixed  attractive/repulsive  schemes  work  better. 

•  finally  anisotropy  has  been  introduced  because  the  previously  defined  scheme  seemed 
well  suited,  and  the  generalization  straightforward. 

A  list  of  desirable  features  has  been  presented  in  the  introduction.  We  now  review  it: 

•  functionality.  The  new  approach  works,  as  the  meshes  in  this  document  testify. 

•  robustness.  Unfortunately  the  scheme  is  not  as  robust  as  it  could  be.  Problems  can 
arise  when  constraints  and  the  feature  size  function  are  “incompatible”,  that  is  if 
the  distance  between  two  constrained  nodes  is  much  shorter  than  the  desired  edge 
length.  The  system  oscillates  as  particles  are  continuously  created  and  annihilated. 
The  system  can  suffer  from  the  same  pathology  when  changes  in  the  metric  tensor  are 
too  abrupt. 

•  quality.  The  quality  of  the  results  is  good,  as  the  angle  and  edge  length  histograms 
show.  The  presented  approach  can  also  be  seen  as  a  mesh  postprocessing  tool.  In  pra¬ 
ctice  it  has  proven  superior  to  Laplacian  smoothing.  The  reasons  seem  to  be  twofold: 
Laplacian  smoothing  doesn’t  have  any  knowledge  of  what  the  inter-point  distance 
should  be,  and  maintenance  of  the  Delaunay  characteristics  after  point  displacement 
improves  the  shape  of  the  mesh  elements. 

•  speed.  It  is  probably  not  a  strength  of  the  presented  approach,  although  the  asym¬ 
ptotic  running  time  seems  to  be  0{n),  where  n  is  the  number  of  nodes  in  the  mesh. 
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However,  in  applications  where  constant  remeshing  is  required  the  new  approach  sho¬ 
uld  perform  well. 

•  minimal  user  interaction.  Several  variables  controlling  the  particle  interaction  model 
need  to  be  set.  Although  an  experienced  user  could  decide  on  what  values  to  use,  it 
is  possible  to  come  up  with  a  standard  set  of  values,  in  which  case  no  user  input  is 
needed  at  this  point. 

•  controllability.  The  introduction  of  a  background  mesh  seems  to  give  the  user  sufficient 
control.  If  he  miscalculates  the  effect  of  the  linear  interpolation,  he  can  always  add  a 
new  sample  to  have  the  mesh  better  fit  his  desires. 


6.1  Future  work 

Directions  for  future  work  are  many: 

•  apply  the  method  to  solve  real  problems.  The  meshes  that  have  been  presented  here 
were  all  generated  with  hand-made  domains  and  hand-made  background  meshes.  The 
addition  of  a  finite  element  solver  would  help  to  demonstrate  the  abilities  of  this  new 
meshing  method  to  generate  adapted  meshes. 

•  generalize  to  non-polygonal  domains.  Most  objects  in  the  real  world  cannot  be  descri¬ 
bed  by  a  set  of  line  segments.  A  more  powerful  scheme,  which  may  include  splines,  is 
thus  needed. 

•  generalize  to  higher  dimensions.  3-dimensional  meshing  is  not  as  well  understood  as 
2-dimensional  meshing.  Maybe  a  physically-based  approach  could  solve  the  sliver 
problem^. 

•  find  even  better  particle  interaction  schemes.  We  believe  that  the  proposed  inte¬ 
raction  scheme  is  good,  but  could  still  be  improved,  especially  on  the  particle  cre¬ 
ation/annihilation  side.  A  scheme  where  particles  are  inserted  at  centers  of  circum- 
circles  could  maybe  generate  better  meshes  that  need  less  smoothing.  On  the  annihi¬ 
lation  side,  collapsing  edges  might  be  better  than  removing  nodes. 

sliver  an  ill-shaped  tetrahedron  of  which  the  four  vertices  are  almost  co-planar 
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•  introduce  an  a  priori  feature  size  function  based  on  geometry  when  none  is  provi¬ 
ded.  This  would  probably  make  the  scheme  more  robust  (see  previous  remark  on 
robustness). 

•  find  out  how  much  remeshing  from  a  previous  solution  is  better  (same  mesh  quality 
at  lower  computational  cost)  than  remeshing  from  scratch. 

6.2  Acknowledgements 

I  would  like  to  thank  Paul  Heckbert  for  inviting  me  to  Carnegie  Mellon,  and  for  being  a 
wonderful  advisor  to  work  with;  Marshall  Bern,  Michael  Garland,  Michael  Erdmann,  Andy 
Witkin,  Jim  Winget,  Scott  Canann,  and  Paul  Heckbert  for  their  helpful  comments  and 
suggestions;  and  my  father  who  has  been  supporting  me  for  more  than  23  years  now. 


55 


56 


Bibliography 


[1]  Eric  B.  Becker,  Graham  F.  Cary,  and  J.  Tinsley  Oden.  Finite  Elements:  An  Introdu- 
etion,  volume  1.  Prentice-Hall,  Englewood  Cliffs,  NJ,  1981. 

[2]  Marshall  Bern  and  David  Eppstein.  Mesh  generation  and  optimal  triangulation.  Techni¬ 
cal  Report  P92-00047,  Xerox  PARC,  1992. 

[3]  M.J.  Castro-Diaz,  F.  Hecht,  and  B.  Mohammadi.  New  progress  in  anisotropic  grid 
adaptation  for  inviscid  and  viscous  flows  simulations.  In  4th  Annual  International 
Meshing  Roundtable,  1995. 

[4]  L.  Paul  Chew.  Constrained  Delaunay  triangulations.  Algorithmiea,  4:97-108,  1989. 

[5]  E.F.  D’Azevedo.  Optimal  triangular  mesh  generation  by  coordinate  transformation.  J. 
Sei.  Stat.  Comput.,  12(4):755-786,  July  1991. 

[6]  L.  Devroye,  E.P.  Miicke,  and  B.  Zhu.  A  note  on  point  location  in  Delaunay  triangula¬ 
tions  of  random  points.  Submited  for  publication,  1995. 

[7]  H.  Edelsbrunner  and  R.  Seidel.  Voronoi  diagrams  and  arrangements.  Dise.  and  Comp. 
Geom.,  8(l):25-44,  1986. 

[8]  D.A.  Field.  Laplacian  smoothing  and  Delaunay  triangulations.  Comm.  Appl.  Num. 
Methods,  4:709-712,  1988. 

[9]  S.  Fortune.  A  sweepline  algorithm  for  Voronoi  diagrams.  Algorithmiea,  2:153-174, 
1987. 

[10]  William  H.  Frey  and  David  A.  Field.  Mesh  relaxation:  A  new  technique  for  improving 
triangulations.  International  Journal  for  Numerieal  Methods  in  Engineering,  31:1121- 
1133,  1991. 

[11]  N.A.  Golias  and  T.D.  Tsiboukis.  An  approach  to  refining  three-dimensional  tetra¬ 
hedral  meshes  based  on  delaunay  transformations.  International  Journal  for  Numerieal 
Methods  in  Engineering,  37:793-812,  1994. 

[12]  James  Gosling  and  Henry  McGilton.  The  Java  Language  Environment:  a  white  paper. 
Sun  Microsystems,  October  1995. 

[13]  Leonidas  Guibas  and  Jorge  Stolfi.  Primitives  for  the  manipulation  of  general  subdi¬ 
visions  and  the  computation  of  Voronoi  diagrams.  ACM  Transaetions  on  Graphies, 
4(2):74-123,  April  1985. 


57 


[14]  L.J.  Guibas,  D.E.  Knuth,  and  M.  Sharir.  Randomized  incremental  construction  of 
Delaunay  and  Voronoi  diagrams.  Technical  Report  STAN-CS-90-1300,  Stanford  Uni¬ 
versity,  1990. 

[15]  K.  Ho-Le.  Finite  element  mesh  generation  methods:  a  review  and  classification. 
Computer-Aided  Design,  20(l):27-38,  Jan/Feb  1988. 

[16]  Thomas  Kao  and  David  M.  Mount.  Dynamic  maintenance  of  Delaunay  triangulations. 
Technical  Report  CS-TR-2585,  University  of  Maryland,  1991. 

[17]  C.  L.  Lawson.  Software  for  C  surface  interpolation.  In  John  R.  Rice,  editor,  Mathe- 
matieal  Software  III,  pages  161-194.  Academic  Press,  1977. 

[18]  Geoff  Leach.  Improving  worst-case  optimal  Delaunay  triangulation  algorithms.  In  fth 
Canadian  Conferenee  on  Computational  Geometry,  1992. 

[19]  D.A.  Lindholm.  Automatic  triangular  mesh  generation  on  surfaces  of  polyhedra.  IEEE 
Trans.  Magneties,  19:2539-2542,  1983. 

[20]  Dani  Lischinski.  Incremental  Delaunay  triangulations.  In  Paul  S.  Heckbert,  editor, 
Graphies  Gems  IV,  pages  47-59.  Academic  Press,  1994. 

[21]  S.H.  Lo.  A  new  mesh  generation  scheme  for  arbitrary  planar  domains.  International 
Jounral  for  Numerieal  Methods  in  Engineering,  21:1403-1426,  1985. 

[22]  Dimitri  J.  Mavriplis.  Adaptive  mesh  generation  for  viscous  flows  using  Delaunay  trian¬ 
gulation.  Journal  of  Computational  Physies,  90(2):271-291,  October  1990. 

[23]  J.  Peraire,  M.  Vahdati,  K.  Morgan,  and  O.C.  Zienkiewicz.  Adaptive  remeshing  for 
compressible  flow  computations.  Journal  of  Computational  Physies,  72:449-466,  1987. 

[24]  Franco  P.  Preparata  and  Michael  Ian  Shamos.  Computational  Geometry:  an  Introdu- 
etion.  Springer-Verlag,  1985. 

[25]  J.  Ruppert.  A  new  and  simple  algorithm  for  quality  2-dimensional  mesh  generation.  In 
4th  ACM-SIAM  Symp.  on  Dise.  Algorithms,  pages  83-92,  1993. 

[26]  Kenji  Shimada.  Physieally-Based  Mesh  Generation:  Automated  Triangulation  of  Sur- 
faees  and  Volumes  via  Bubble  Paeking.  PhD  thesis,  MIT,  1993. 

[27]  Kenji  Shimada  and  David  C.  Gossard.  Computational  methods  for  physically-based  fe 
mesh  generation.  In  PROLAMAT.  IFIP  TC5/WG5.3,  1992. 

[28]  I.S.  Sokolnikoff.  Tensor  Analysis,  Theory  and  Applieations  to  Geometry  and  Meehanies 
of  Continua.  John  Wiley,  New  York,  2nd  edition,  1964. 

[29]  Richard  Szeliski  and  David  Tonnesen.  Surface  modeling  with  oriented  particle  systems. 
In  SIGGRAPH’92  Proeeedings,  pages  185-194,  1992. 

[30]  Greg  Turk.  Generating  textures  on  arbitrary  surfaces  using  reaction-diffusion.  In 
SIGGRAPH’Ol  Proeeedings,  pages  289-298,  1991. 

[31]  Marie-Gabrielle  Vallet.  Generation  de  maillages  anisotropes  adaptes  -  application  a  la 
capture  de  couches  limites.  Technical  Report  1360,  INRIA,  1990. 


58 


[32]  N.P.  Weatherill  and  O.  Hassan.  Efficient  three-dimensional  Delaunay  triangulation  with 
automatic  point  creation  and  imposed  boundary  constraints.  International  Journal  for 
Numerical  Methods  in  Engineering,  37:2005-2039,  1994. 

[33]  William  Welch.  Serious  Putty:  Topological  Design  for  Variational  Curves  and  Surfaces. 
PhD  thesis,  Carnegie  Mellon  University,  1995. 

[34]  Andrew  P.  Within  and  Paul  S.  Heckbert.  Using  particles  to  sample  and  control  implicit 
surfaces.  In  SIGGRAPH’Of  Proceedings,  1994. 

[35]  M.A.  Yerry  and  M.S.  Shepard.  A  modified  quadtree  approach  to  finite  element  mesh 
generation.  IEEE  Computer  Graphics  and  Applications,  3:39-46,  January /February 
1983. 


59 


